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4.  INTRODUCTION 


This  project  is  to  facilitate  research  in  digital  mammography  and  related  technologies,  in 
particular  computer-aided  diagnosis  (CAD)  and  image  processing.  A  major  limitation  to  the 
rapid  development  and  subsequent  clinical  implementation  of  these  technologies  is  the  lack  of  a 
standardized  set  of  mammograms  with  absolute  truth  (e.g.,  exact  location  and  extend  of  a  cancer) 
to  be  used  in  development  and  evaluation.  We  are  developing  a  method  to  produce  computer- 
simulated  mammograms.  The  approach  is  to  model  the  creation  of  the  mammogram  on  the 
computer  —  all  steps  from  x  rays  exiting  the  breast  to  the  image  being  displayed  on  a  light  box. 
The  basic  model,  which  we  have  developed  previously,  has  been  improved  and  is  combined  with 
accurate  information  of  the  appearance  of  normal  breast  anatomy  and  of  benign  and  malignant 
breast  lesions.  These  are  obtained  from  high  quality  images  of  cadaver  breasts  and  biopsy 
specimens.  Our  approach  is  similar  to  that  of  Van  Metter  et  al.  who  modeled  chest  radiographs 
[Van  Metter  1986].  We  believe  that  this  technique  can  produce  simulated  mammograms  that 
appear  to  be  actual  mammograms.  This  hypothesis  is  being  tested  by  performing  quantitative 
comparisons  of  simulated  and  real  mammograms. 

5.  BODY  OF  REPORT 
5.1  Tasks 

We  originally  proposed  four  tasks  in  the  Statement  of  Work,  which  are  listed  below.  Not 
all  the  tasks  were  completed  and  some  of  the  tasks  were  modified. 

1 .  To  obtain  radiographs  of  mastectomy  and  tissue  specimens: 

(a)  radiograph  100  different  mastectomy  breast  tissues  at  2.0  times  geometric  magnification 
recording  image  on  direct  film  (without  intensifying  screen)  at  five  different  orientations; 

(b)  radiograph  240  different  tissue  specimens  at  4.0  times  geometric  magnification  recording 
images  on  direct  film  (without  intensifying  screen)  at  five  different  orientations;  and 

(c)  segment  lesions  from  specimen  radiographs  and  measure  their  size,  contrast,  and  shape 
metrics. 

2.  To  develop  further  a  computer  model  of  image  formation: 

(a)  modify  previously  developed  model  for  point  source  versus  parallel  beam; 

(b)  measure  and  model  detector  noise  for  film  digitizer  and  screen-film  system; 

(c)  measure  scatter  as  a  function  of  position  in  the  image;  and 

(d)  measure  beam  intensity  as  a  function  of  position  in  the  image. 

3.  To  produce  simulated  mammograms: 

(a)  produce  simulated  mammograms  with  and  without  lesions; 
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(b)  make  preliminary  comparison  to  actual  mammograms;  and 

(c)  make  adjustments  to  model,  if  necessary. 

4.  To  evaluate  simulated  images: 

(a)  collect  real  mammograms:  normals  and  those  with  lesions; 

(b)  compare  real  and  simulated  mammograms  based  on  quantitative  measurements; 

(c)  conduct  pilot  observer  study;  and 

(d)  conduct  observer  study  comparing  ROIs  from  real  and  simulated  mammograms. 


5.1.1  Obtain  radiographs  of  tissue  specimens  and  mastectomies 

We  partially  completed  tasks  (a)  and  (c),  but  we  did  not  start  (b).  We  have  decided  to  use 
cadaver  breasts  instead  of  mastectomy  specimens,  because  they  are  easier  to  obtain.  We  have 
fourteen  cadaver  breasts  that  have  been  imaged  using  plain  film.  We  have  not  yet  collected 
specimen  radiographs.  Figure  1  shows  an  example  of  a  piece  of  a  cadaver  breast. 

5.1.2  Further  development  of  computer  simulation  method 

The  theoretical  basis  for  the  model  has  been  developed  previously  by  the  PI,  but  with  a 
number  of  simplifying  assumptions  [Nishikawa  1989,  1990a,  b].  For  this  project,  we  need  to 
check  these  assumptions  and  include  other  relevant  factors  particular  to  our  application.  In 


Figure  1 .  A  portion  of  a  screen-film  image  of  a  cadaver  breast.  Cadaver  breasts  are  imaged  with  high 
fidelity  (film  only,  no  screen,  using  a  high  x-ray  exposure  and  geometric  magnification).  Those  images 
are  digitized  at  50  microns  and  used  as  input  to  our  simulation  model.  Folds  in  the  skin  are  apparent  in 
the  image  and  can  be  removed  by  placing  the  breast  in  saline. 
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addition,  we  need  to  covert  the  theory  into  a  computer  program  that  can  produce  a  simulated 
image.  Our  efforts  in  these  areas  are  described  below. 

We  originally  planned  to  modify  our  model  to  account  for  scattered  radiation  and  non¬ 
uniformity  of  the  x-ray  beam.  We  still  plan  to  incorporate  these  modifications,  as  they  are  fairly 
straightforward  to  implement.  However,  instead  of  working  on  these  modifications,  we  chose 
instead  to  incorporate  k-fluorescence  re-absorption  into  our  model  as  this  is  not  straightforward 
and  greatly  expands  the  usefulness  of  our  model.  At  the  time  of  submission  of  our  proposal, 
Gd02S2  was  the  dominant  phosphor  in  use  for  mammography,  both  screen-film  systems  and 
experimental  (full-field)  and  commercial  (small-field)  digital  systems.  The  k-edge  of  this 
phosphor  is  well  above  the  energies  used  in  mammography,  so  k-fluorescence  will  not  occur. 
However,  the  dominant  phosphors  in  digital  systems  today  are  Csl  and  selenium.  Csl  has  a  k- 
edge  at  approximately  33  and  36  keV  and  selenium  at  approximately  13  keV.  These  energies  are 
important  for  digital  mammography,  where  x-ray  beams  as  high  as  40  kVp  are  used  in 
commercial  systems.  Therefore,  to  be  accurate,  k-fluorescence  re-absorption  needs  to  be 
incorporated  into  our  model. 

5.1.2.a.  Modify  model  from  parallel  beam  of  x  rays  to  x  rays  from  a  point  source. 

In  our  previous  model,  we  assumed  that  the  x  rays  incident  on  the  detector  were  at  normal 
incidence  (parallel  beam).  In  practice,  most  x  rays  are  incident  at  angles  slightly  less  than  90 
degrees.  As  illustrated  in  Fig.  2,  x  rays  originating  from  a  focal  spot  a  finite  distance  from  the 
detector  will  produce  an  x-ray  beam  that  is  diverging  and  is  incident  on  the  detector  at  some 
angle  defined  by  0,  where  0  depends  on  the  distance  of  the  focal  spot  from  the  detector,  and  the 
x-y  position  of  the  x  ray  incident  on  the  detector  -  the  larger  the  x-y  position,  the  larger  the  angle. 
If  the  x  rays  are  parallel  and  are  incident  at  90  degrees  to  the  detector,  then  each  pixel  at  the 
output  has  a  one-to-one  correspondence  to  a  point  at  the  input,  in  terms  of  energy  deposition 
(ignoring  scattered  x  rays).  If  the  beam  diverges  however,  depending  on  the  depth  of  interaction 
of  the  x  ray  in  the  detector,  a  pixel  at  the  output  can  have  contributions  from  x  rays  that  entered 
the  detector  at  several  different  points  at  the  input  surface.  In  other  words,  the  position  of  where 
the  x  ray  would  actually  interact  will  be  different  from  the  position  calculated  by  the  model  and 
this  difference  is  A. 

We  have  calculated  the  difference  in  energy  deposition  in  the  detector  assuming  a  parallel 
beam  and  a  diverging  x-ray  beam.  Standard  patient  geometry  is  a  focal-film  distance  of  65  cm 
and  maximum  detector  size  of  24x30  cm,  with  the  central  axis  of  the  beam  at  one  edge  of  the 
detector,  centered  in  the  other  direction).  Then  assuming  a  phosphor  thickness  of  85  microns,  we 
have  calculated  that  the  maximum  energy  spread  because  of  a  diverging  beam  is  35  microns  (i.e., 
the  difference  between  parallel  and  non-parallel  beam  assumptions  is  a  35-micron  difference  in 
the  point  at  which  the  x  ray  is  absorbed).  This  will  occur  at  the  periphery  of  the  detector  and  only 
for  those  x  rays  interacting  in  the  bottom  of  the  phosphor  (i.e.,  furthest  from  the  source). 
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Figure  2.  A  schematic  representation  of  the  geometry  of  the  x-beam  incident  on  the  phosphor  screen  as 
seen  from  a  coronal  view.  The  x-ray  beam  originates  from  a  focal  spot  that  is  65  cm  above  the  screen. 
As  a  result,  most  x  rays  are  incident  on  the  detector  at  some  angle  0>O.  In  our  model,  we  assume  that 
the  x  rays  are  incident  at  0=0  degrees  (normal  incidence)  as  if  they  were  a  parallel  beam  of  x-rays.  This 
can  produce  a  discrepancy  (A)  between  where  the  x  ray  actually  interacts  in  the  screen  and  where  it 
interacts  according  to  our  model.  The  value  of  A  depends  on  the  depth  at  which  the  x  ray  interacts  and 
0,  which  depends  on  the  x-y  position  of  the  incident  x  ray.  Note  that  the  compression  plate  is  not 
shown  and  that  the  angle  depicted  has  been  exaggerated  to  illustrate  the  point. 


Phosphor  Screen 


X-ray 
Beam 


Figure  3  shows  the  fraction  of  x  rays  entering  the  phosphor  at  a  given  x-y  location  that 
will  interact  in  the  screen  more  than  17  microns  laterally  from  the  point  it  entered  the  screen. 
That  is,  this  figure  shows  the  fraction  of  x-ray  interactions  that  will  give  rise  to  more  than  a  17- 
micron  error  when  assuming  all  x-rays  enter  the  screen  at  right  angles  compared  to  the  real 
situation  where  the  x  rays  enter  the  screen  at  some  angle  slightly  less  than  90  degrees.  For  this 
calculation,  it  is  assumed  that  the  x  rays  originate  from  a  point  source  65  cm  above  the  page, 
centered  at  0,0.  For  example,  a  17-micron  difference  will  occur  at  the  very  farthest  comer  of  the 
screen  for  less  than  40%  of  all  x-ray  interactions  at  that  point.  The  average  size  compressed 
breast  is  approximately  100  cm2.  Therefore,  for  most  breasts  the  difference  between  a  parallel 
beam  of  x  rays  and  a  diverging  beam  is  negligible.  For  larger  breasts,  the  difference  is  minimal. 

Since  we  will  initially  use  a  17-micron  pixel  size  at  the  output  to  form  the  image,  and 
subsequently  form  larger  pixels  by  averaging,  we  need  to  take  the  beam  divergence  into  account 
to  be  absolutely  accurate.  However,  only  12.5%  of  all  interactions  in  the  screen  will  lead  to  a 
small  error  and  most  if  not  all  of  these  will  be  for  pixels  outside  of  the  breast  area.  Because  of 
the  added  complexity  a  non-parallel  beam  introduces,  we  have  initially  assumed  a  parallel  beam 
in  our  calculations.  If  we  find  that  this  assumption  produces  poor  results,  we  will  modify  our 
calculations  for  a  diverging  beam. 
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Figure  3.  This  iso-plot  shows  the  fraction  of  x-rays  that  interact  within  a  phosphor  screen  for  which  there 
is  a  17-micron  difference  between  two  x-rays  beams  entering  the  phosphor  at  the  same  location  one  at 
normal  incidence  and  one  at  angle.  The  x  rays  are  assumed  to  have  originated  from  a  point  65  cm  above 
the  surface  of  the  screen  directly  above  the  (0,0)  point  of  the  screen  and  the  screen  is  85-microns  thick. 

The  typical  screen-film  system  used  in  mammography  is  1 8x24  cm.  For  larger  breasts,  a  24x30  cm  system 
is  used.  The  inner  contour  covers  an  area  of  265  cm2,  which  is  more  than  twice  the  size  of  an  average 
compressed  breast. 


Note  that  the  image  is  still  formed  using  a  diverging  beam,  so  that  the  appearance  of  the 
tissue  is  correct.  It  is  only  once  that  the  x  rays  impinge  on  the  screen  that  we  assume  that  they 
enter  the  screen  at  normal  incidence,  even  though  they  may  have  arrived  at  the  screen  at  not 
normal  incidence. 

5.1.2.b.  Model  detector  noise  for  film  digitizer  and  screen-film  system. 

Based  on  published  noise  power  spectra  of  film  [Bunch  1 999]  we  integrated  the  spectra, 
weighted  by  the  Fourier  spectrum  of  a  50-micron  scanning  aperture.  This  gives  us  the  standard 
deviation  (square  root  of  the  integral)  as  a  function  of  film  density  (see  Figure  4).  This  was  done 
for  different  film  optical  densities.  Then  assuming  the  noise  on  the  film  follows  Gaussian 
statistics,  we  used  a  random  number  generator  to  produce  a  noise  pattern  that  is  added  to  the 
simulated  image.  That  is,  given  a  pixel  with  a  certain  film  density,  we  generate  a  random  number 
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Figure  4.  The  relationship  between  the  noise  for  a  50-micron  pixel  (as  given  by  the  standard  deviation) 
as  a  function  of  the  film  optical  density.  This  gives  the  magnitude  of  the  noise  due  to  film  granularity. 


that  is  from  a  Gaussian  distribution  with  zero  mean  and  a  standard  deviation  corresponding  to 
that  film  density,  and  then  add  this  number  to  the  pixel  value. 

We  have  also  data  on  the  noise  of  the  film  digitizer  as  a  function  of  film  density  that  we 
have  measured  in  our  laboratory  (see  Figure  5).  To  obtain  these  data,  we  digitized  a  calibration 
strip  that  had  near-noiseless  squares  of  constant  optical  densities  from  0.10  to  4.0.  Again, 
assuming  Gaussian  statistics,  a  random  number  generator  is  used  to  produce  a  noise  pattern  that 
is  added  to  the  simulated  image.  This  information  is  necessary  because  we  are  using  digital  data 
for  our  comparison  (simulated  versus  real)  observer  study,  so  we  need  to  model  digitization 
noise.  We  are  using  digital  data  because,  if  we  printed  the  simulated  images  on  film  and  used 
film  for  comparison,  it  would  be  obvious  which  images  were  simulated  because  they  would  look 
pixilated.  Note  that  compare  to  film  granularity  noise,  the  noise  due  to  the  digitizer  is 
insignificant  under  most  conditions.  However,  at  film  densities  above  approximately  3.0,  the 
two  noise  sources  are  comparable. 

5.1 .2.c.  Measure  scatter  as  a  function  of  position  in  the  image 

Our  original  plan  was  to  measure  the  scatter  as  a  function  of  position  in  the  image  by 
placing  lead  blockers  under  the  excised  breast  specimen.  This  would  allow  us  to  calculate  the 
scatter  to  primary  ratio  at  several  points  in  the  breast.  Then  by  interpolation,  we  could  estimate 
the  scatter  field  for  the  whole  breast.  This  scatter  field  would  then  be  added  to  the  simulated 
image.  The  drawback  of  this  method  is  that  we  would  need  to  repeat  the  experiment  for  every 
breast  specimen  that  we  would  be  simulating  and  the  result  would  not  necessarily  be  accurate, 
because  of  uncertainties  in  interpolation. 
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Figure  5.  The  noise  power  spectrum  of  the  film  digitizer  for  a  film  with  an  optical  density  of  1 .5. 


Instead,  we  will  incorporate  scatter  into  our  model  using  a  scatter  term  based  on  measured 
point-spread  functions  for  scatter  measured  by  Boone  et  al.  [Boone  2002].  For  each  pixel  at  the 
output  of  the  detector,  the  number  of  x  rays  interacting  in  a  column  through  the  detector  is 
computed.  That  is,  at  the  output  of  the  phosphor,  there  is  an  image  corresponding  to  the  number 
of  x  rays  interacting  in  the  phosphor  at  each  pixel.  We  will  convolve  this  with  the  point-spread 
function  for  scatter,  which  is  dependent  on  the  breast  thickness  and  very  weakly  on  the  energy  of 
the  x-ray  beam.  Figure  6  shows  a  point-spread  function  of  scatter  for  a  30-kVp  beam  and  a  6-cm 
thick  breast. 

We  have  not  yet  incorporated  scatter  into  our  model,  as  we  are  working  on  the  noise 
aspects  of  the  model  first.  It  will.be  straightforward  to  implement  because  we  only  need  to 
perform  a  convolution.  This  will  be  done  in  the  spatial  frequency  domain  as  a  multiplication  by 
taking  the  Fourier  transform  of  the  scatter  point-spread  function  and  the  image  at  the  output  of 
the  detector.  We  already  perform  some  of  the  calculations  of  the  model  in  the  spatial  frequency 
domain,  so  the  image  is  already  in  the  spatial  frequency  domain.  We  do  not  anticipate  any 
difficulty  implementing  this  correction. 

5.1 .2.d.  Measure  beam  intensity  as  a  function  of  position  in  the  image 

The  intensity  distribution  of  x  rays  from  the  Faxitron  does  not  match  that  of  a  standard  x- 
ray  machine.  Therefore,  to  more  accurately  simulate  mammograms  this  difference  needs  to  be 
corrected.  This  is  done  by  multiplying  the  input  image  to  our  simulation  program  by  the  ratio  of 
the  x-ray  intensity  of  a  mammography  machine  to  the  x-ray  intensity  of  the  Faxitron.  We  have 
measured  the  x-ray  intensity  of  the  Faxitron,  but  we  haven’t  measure  the  intensity  on  a 
mammography  system,  yet.  Since  we  are  only  simulating  limited  sized  areas  within  a 
mammogram,  this  correction  is  a  second  order  one.  Further,  at  this  stage  we  are  only  comparing 
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Figure  6.  The  point-spread  function  due  to  scattered  radiation  within  the  patient  for  compressed  breast 
thicknesses  ranging  from  2  cm  to  8  cm.  The  inset  shows  the  2D  distribution  of  the  point  spread 
function  for  a  6-cm  compressed  breast.  The  data  are  for  a  30  kVp  mammographic  spectrum,  although 
the  point  spread  functions  are  nearly  independent  of  kVp. 

images  taken  on  the  Faxitron.  We  do  not  anticipate  any  difficulty  implementing  this  correction 
in  the  future. 

5.1 ,2.e.  K-fluorescence  re-absorption  (research  beyond  original  statement  of  work) 

Many  x-ray  detectors  use  a  phosphor  screen  to  convert  the  x-ray  quanta  into  visible  light 
photons  that  are  subsequently  collected  by  an  optical  detector.  Fluorescent  x  rays  can  be 
produced  within  the  phosphor.  Some  of  these  fluorescent  x-ray  quanta  are  reabsorbed  in  the 
phosphor.  This  fluorescent  x-ray  re-absorption  affects  the  quality  of  images.  We  have  developed 
a  3D  convolution  model  to  numerically  simulate  the  characteristic  x-ray  re-absorption  inside  an 
image  detector  or  intensifying  screen  with  Monte  Carlo  technology.  The  3D  model  is  built 
through  analyzing  the  physic  process  of  creation,  traveling,  and  absorption  of  characteristic  x-ray 
quanta  so  that  it  can  simulate  the  noise  properties  of  characteristic  x-ray  re-absorption  accurately. 
The  phosphor  layer  of  an  image  detector  is  conceptually  divided  into  multiple  sub-layers  and 
every  sub-layer  is  conceptually  divided  into  multiple  voxels  based  on  the  pixel  size  of  a  digital 
image.  The  3D  model  computes  the  number  distribution  of  reabsorbed  characteristic  x-ray 
quanta  in  every  sub-layer  using  a  pre-calculated  convolution  kernel.  The  convolution  kernel 
gives  the  probability  distribution  that  a  characteristic  x-ray  quantum  can  be  reabsorbed  by  a 
voxel.  A  cut-off  window  of  the  convolution  is  calculated  to  reduce  the  computation  time.  This 
3D  convolution  model  shows  that  most  characteristic  x-ray  re-absorption  blurs  the  image. 
However,  some  of  characteristic  x-ray  re-absorption  enhances  the  primary  x-ray  input,  depending 
on  the  relative  position  between  the  creation  site  and  re-absorption  site  of  the  fluorescent  x  rays. 
The  simulation  result  is  used  to  quantitatively  analyze  the  effect  of  the  fluorescent  x-ray  re- 
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absorption  on  the  imaging  performance  of  the  detector  or  intensifying  screen.  The  calculation  to 
a  cesium  iodide  layer  with  pixel  size  17x17  microns  shows  that  the  characteristic  re-absorption 
may  reduce  the  modulation  transfer  function  of  the  image  detector  as  much  as  10%  at  some 
spatial  frequency,  which  is  close  to  the  measurement.  It  also  reduces  the  NPS  at  non-zero  spatial 
frequencies. 

The  details  of  the  method  are  given  in  a  preprint,  which  is  attached  to  this  report.  We 
also  develop  one  new  scientific  finding  in  this  work.  We  have  demonstrated  that  the  depth 
within  the  phosphor  at  which  the  fluorescent  x-ray  is  re-absorbed  affects  image  quality.  Metz 
and  Vybomy  previously  showed  theoretically  that  the  effect  on  the  NPS  of  the  re-absorption  is  to 
decrease  the  NPS  at  high  spatial  frequencies  by  a  constant  amount  [Metz  1983].  This  work  has 
been  confirmed,  again  theoretically  using  different  methods,  by  Hillen  et  al.  [Hillen  1991]  and 
Yao  and  Cunningham  [Yao  2001].  In  the  theoretical  models  used  to  date,  it  is  assumed  that  the 
phosphor  screen  is  one  slab  of  phosphor.  In  our  work,  we  assume  the  phosphor  is  made  of  a 
stack  of  many  thin  layers  of  phosphor.  The  affect  this  has  on  the  NPS  is  shown  in  Fig.  7.  The 
results  show  that  when  the  screen  is  properly  modeled  the  reduction  in  the  NPS  at  high  spatial 
frequency  is  much  more  dramatic.  In  the  single  slab  model,  the  reduction  is  20%.  When  the 
screen  is  modeled  as  many  thin  slabs,  our  results  show  that  the  reduction  is  70%. 


Figure  7.  This  shows  the  ratio  of  the  shape  of  the  noise  power  spectmm  (NTF)  when  the  effect  of  K- 
fluorescence  is  modeled  to  when  it  is  not.  The  different  curves  show  the  effect  of  the  depth  within  the 
phosphor  that  the  K-fluorescence  is  reabsorbed.  As  the  screen  is  considered  to  be  composed  of  more 
layers,  the  ratio  decreases  at  high  spatial  frequencies.  Previous  work  in  this  area  considered  the  screen 
to  be  a  single  thick  slab,  resulting  in  an  underestimate  of  the  effect  of  k-fluorescence  reabsorption  on 
the  noise  power  spectmm.  The  curves  were  generated  by  using  a  Monte  Carlo  model. 
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5.1.3.  Produce  simulated  mammogram 

We  have  tested  our  simulation  technique  using  phantoms  and  standard  image  quality 
metrics:  modulation  transfer  function  (MTF)  and  noise  power  spectrum  (NPS).  Since  images 
will  be  made  on  the  Kodak  Min-R  2000  screen-film  system,  we  have  compared  the 
characteristics  of  our  simulated  images  with  those  for  the  Min-R  2000  system,  based  on 
published  data. 

To  examine  the  resolution  properties  of  our  simulated  image,  we  simulated  a  point  source 
of  x  rays  (see  Fig.  8)  and  used  it  to  calculate  the  modulation  transfer  function  (MTF)  (see  Fig.  9). 
The  MTF  from  the  simulated  image  and  published  MTF  [Bunch  1 999]  show  excellent 
agreement  -  better  than  5%. 

To  examine  the  noise  properties  of  our  simulated  image,  we  generated  a  uniform  random 
noise  image  using  a  Gaussian  random  number  generator  (see  Figure  10).  The  output  image  from 
our  simulation  model  is  then  used  to  calculate  the  NPS.  Currently,  we  do  not  have  accurate 
information  of  the  conversion  from  number  of  light  quanta  incident  on  the  film  to  optical  density. 
This  will  affect  the  magnitude  of  the  NPS,  but  not  its  shape.  Therefore,  we  compared  the  shapes 
of  the  NPS  for  our  simulated  image  and  for  the  Min-R  2000  system  (see  Figure  11).  There  is 
good  agreement  -  less  than  8%  difference  except  at  high  spatial  frequencies.  Since  we  have  not 
included  film  granularity  in  our  simulated  image,  we  would  expect  that  our  simulated  image 
would  have  lower  NPS  at  high  spatial  frequencies. 

We  have  also  made  simulated  images  of  the  phantoms.  The  phantoms  were  first  imaged 
on  Kodak  XV  film  (non-screen  system)  using  three  times  geometric  magnification  that  will 
produce  a  high  fidelity  image,  with  low  noise.  The  images  taken  on  the  XV  film  are  digitized  at 
50-micron  pixel  resolution.  This  produces  an  image  with  effective  pixel  size  of  1 7  microns.  The 
resulting  image  is  put  through  our  simulation  model  to  produce  an  image  with  the  same  exposure 
as  used  for  the  standard  mammography  system  -  a  Kodak  Min-R  2000  system. 

To  compare  the  simulated  image  to  the  screen-film  image,  we  digitized  the  screen-film 
image  at  50  microns  and  construct  a  50-micron  resolution  simulated  image  (see  Figure  12). 
Visually  the  real  and  the  simulated  images  look  similar,  except  that  the  simulated  image  is 
slightly  darker  than  the  real  image.  This  is  caused  by  using  an  incorrect  conversion  from  number 
of  light  quanta  to  film  optical  density.  We  are  now  in  the  process  of  correcting  this  problem  after 
discussion  with  Dr.  Phillip  Bunch  of  Eastman  Kodak  Company. 
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Figure  8.  The  input  point  image  is  shown  on  the  left,  the  resulting  output  of  our  simulation  model 
shown  on  the  right.  This  point  spread  function  is  then  subjected  to  a  2D  fast  Fourier  transform  to 
calculate  the  modulation  transfer  function  (MTF)  for  the  modeled  phosphor  screen  (Min-R  2000). 


Spatial  Frequency  (cycles/mm) 


Figure  9.  MTF  curves  for  a  Kodak  Min-R  2000  screen-film  system  (x’s),  from  published  data  [Bunch 
1999]  and  from  the  simulated  imaged  (solid  line)  shown  in  Fig.  8. 
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Figure  10.  The  image  on  the  left  contains  uniform  uncorrelated  Gaussian  noise  and  is  used  as  input 
to  our  simulation  model.  The  image  on  the  right  is  the  output  of  the  model.  This  image  and  others 
like  it  are  used  to  calculate  the  noise  power  spectrum. 


Figure  1 1 .  The  image  on  the  left  is  the  2D  noise  power  spectrum  of  the  output  of  the  simulation 
model  for  a  uniform  noisy  input  (see  Figure  10).  The  graph  shows  a  slice  through  the  2D  noise 
power  spectrum  (black  triangles)  and  compares  the  simulated  noise  power  spectrum  to  a  Kodak  Min- 
R  2000  system  provided  to  us  by  P.  Bunch  (red  x’s).  There  is  good  agreement  up  to  approximately  5 
cycles/mm.  The  discrepancy  at  higher  spatial  frequencies  is  because  the  simulation  result  does  not 
include  the  effects  of  film  noise  as  the  experimental  result  does. 
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Figure  12.  Comparison  of  real  image  (left  side)  and  simulated  phantom  image  (right  side)  of  a  portion 
of  a  contrast-detail  phantom.  The  real  image  was  made  using  a  Kodak  Min-R  2000  screen-film  system 
and  was  digitized  at  50-micron  pixel  resolution 


5.1.4.  Evaluation  of  simulated  mammograms 

This  has  not  yet  been  accomplished,  but  we  plan  to  test  whether  radiologists  can  tell  the 
difference  between  our  simulated  and  real  mammograms.  We  will  perform  an  observer  study 
involving  8  radiologists  in  a  two-alternative  forced  choice  (2-AFC)  experiment  [Green  1974]. 
The  radiologists  will  view  a  pair  of  mammograms,  actually  sections  of  a  mammogram.  One  will 
be  a  simulated  image  and  one  will  be  a  real  image.  The  radiologist  will  have  to  choose  which 
one  is  the  real  mammogram.  To  determine  the  number  of  cases  needed,  we  will  first  perform  a 
pilot  study  using  100  cases  and  3  medical  physicists  who  have  worked  in  mammography  for 
more  than  7  years.  We  will  use  standard  methodology  for  performing  and  analyzing  the 
experiment  [Green  1974;  Loo  1984]. 

The  images  will  be  viewed  on  a  pair  of  calibrated  and  identical  CRT  monitors.  It  would 
be  ideal  to  compare  the  original  screen-film  mammogram  to  the  simulated  image  printed  on  film. 
However,  we  believe  that  the  printed  image  would  look  obvious  compared  to  the  real 
mammogram  and  thus  bias  the  experiment. 


5.2.  Recommendations  in  relation  to  the  Statement  of  Work 

There  were  two  changes  to  our  statement  of  work.  The  first  was  to  use  cadaver  breasts 
instead  of  mastectomy  samples,  since  it  is  easier  to  get  cadaver  breasts.  The  second  is  we 
investigated  the  role  of  k- fluorescence  reabsorption  in  our  model.  This  physical  phenomenon  is 
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important  in  the  newer  phosphors  used  in  digital  mammography. 


6.  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Method  for  simulating  mammograms  developed 

•  Linear  systems  model  of  a  phosphor-based  detector  has  been  refined. 

•  New  understanding  of  the  effect  of  k-fluorescence  reabsorption  on  the  noise  power  spectrum 
of  phosphor-based  detectors 


7.  REPORTABLE  OUTCOMES 

We  have  presented  parts  of  this  work  at  one  national  and  three  international  conferences.  We 
have  published  one  conference  proceeding  (attached)  and  we  are  working  on  three  papers  for 
peer-reviewed  journals  (2  are  attached). 

Oral  Presentation 

1.  Nishikawa  RM  and  Lee  S:  A  Method  for  Producing  Simulated  Mammograms.  Presented  at 
the  Presented  at  the  42nd  Annual  Meeting  of  the  American  Association  of  Physicists  in 
Medicine,  July  2002,  Montreal  PQ. 

2.  Nishikawa  RM  and  Lee  S:  A  method  for  producing  simulated  mammograms.  Presented  at  the 
7th  International  Conference  on  Digital  Mammography,  June  2002,  Bremen,  Germany 

3.  Zhang  Y,  Nishikawa  RM:  Computer  Simulation  of  Mammographic  Imaging  for  Applications 
in  CAD.  Presented  at  5th  International  Workshop  on  Computer-Aided  Diagnosis.  June  2004, 
Chicago,  IL. 

4.  Chinander  MA,  Nishikawa  RM,  Zhang  Y:  Numerical  simulation  of  mammographic  imaging 
system  and  assurance  checking.  Presented  at  90th  Scientific  Assembly  and  Annual  Meeting  of 
the  Radiological  Society  of  North  America,  November  2004,  Chicago,  IL. 

Conference  Proceedings 

1.  Zhang  Y,  Nishikawa  RM:  Computer  Simulation  of  Mammographic  Imaging  for  Applications 
in  CAD.  Proc.  CARS  2004  (in  press). 

Manuscripts  in  Preparation  for  peer-reviewed  Journal 

We  are  preparing  3  manuscripts  to  be  submitted  for  peer-review.  Two  are  in  intermediate  draft 
form  and  are  attached.  We  are  still  working  on  the  first  draft  of  the  third  paper. 
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1.  Zhang  Y,  Nishikawa  RM:  Numerical  Simulation  of  Mammographic  Imaging  System.  To  be 
submitted  to  Medical  Physics,  (draft  attached) 

2.  Zhang  Y,  Nishikawa  RM:  3D  Modeling  and  Simulation  of  Characteristic  x-ray  re-absorption 
by  Monte  Carlo  Simulation.  To  be  submitted  to  Computers  in  Medicine  and  Biology,  (draft 
attached) 

3.  Nishikawa  RM,  Zhang  Y,  Chinander  MA:  Effect  of  the  finite  thickness  of  a  phosphor  screen 
and  k-fluorescence  reabsorption  on  image  quality.  To  be  submitted  to  Medical  Physics,  (in 
preparation) 


8.  CONCLUSIONS 

We  are  now  capable  of  producing  simulated  mammograms.  The  resolution  and  noise 
properties  are  closely  matched  to  those  of  real  images.  We  are  still  improving  the  model  with  the 
addition  of  film  noise,  scattered  radiation,  and  light  quanta  to  film  optical  density  conversion.  In 
addition,  we  have  included  the  effects  of  k-fluorescence  reabsorption  in  our  model.  This  will 
allow  us  to  better  simulate  a  wider  variety  of  detectors  used  in  mammography. 

We  have  also  shown  for  the  first  time  that  the  finite  thickness  of  the  screen  will  increase 
the  effect  of  k-fluorescent  re-absorption.  When  the  phosphor  is  treated  as  a  single  thick  slab,  k- 
fluorescence  re-absorption  reduces  the  noise  power  spectrum  by  about  20%.  When  the  phosphor 
is  more  accurately  modeled  as  being  a  stack  of  many  thin  layers,  k-fluorescence  re-absorption 
reduces  the  noise  power  spectrum  by  about  70%. 

As  was  in  our  original  statement  of  work,  we  still  plan  to  conduct  an  observer  study  to 
show  that  the  simulate  images  look  like  real  mammographic  images.  Currently,  we  have  funding 
through  a  training  grant  for  a  post-doctoral  fellow  to  continue  the  work  for  one  year.  We  will 
conduct  the  observer  study  and  make  further  refinements  to  the  model. 
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Numerical  Simulation  of  Mammographic  Imaging  System 

Yinghui  Zhang  and  Robert  M.  Nishikawa 
Department  of  Radiology,  the  University  of  Chicago 
Abstract:  A  method  is  described  for  the  computer  simulation  of  mammographic  imaging. 
Software  has  been  developed  to  simulate  the  physics  processes  of  converting  input  x-ray 
quanta  to  visible  light  photons  inside  a  fluorescent  intensifying  screen.  The  simulation 
algorithm  is  based  on  the  model  of  the  spatial-frequency-dependent  detective  quantum 
efficiency  of  fluorescent  screens,  which  was  developed  previously.  The  x-ray  input  to  the 
simulation  software  was  derived  from  images  taken  on  directly  contacting  film  by 
calibrating  the  film  optical  density  to  the  x-ray  exposure  incident  to  the  film.  The  physics 
processes  that  convert  x-ray  quanta  to  visible  light  photons  inside  a  fluorescent  screen 
were  simulated  by  five  mathematical  transformations.  Standard  FFT  techniques  were 
used.  Swank’s  one-dimensional  fluorescent  screen  modulation  transfer  function  model 
was  extended  into  a  two-dimensional  photon  diffusion  model  by  solving  a  2-dimensional 
Boltzmann  diffusion  equation.  During  simulation  computation,  the  intensifying 
fluorescent  screen  was  conceptually  divided  into  multiple  thinner  sub-layers.  The  final 
image  was  obtained  by  integrating  the  visible  light  output  of  every  sub-layer  over  the 
whole  intensifying  screen.  The  software  based  on  the  simulation  algorithm  in  this  paper 
was  tested  for  assurance  by  three  comparisons:  phantom  images,  modulation  transfer 
function,  and  noise  power  spectrum.  All  test  results  showed  that  the  simulation  method 
and  the  simulation  software  worked  correctly.  Qualitatively,  the  simulated  and  actual 
images  of  test  objects  looked  similar.  For  the  modulation  transfer  function  and  the  shape 
of  the  noise  power  spectra,  the  difference  between  the  simulation  and  the  measured 
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values  of  an  actual  fluorescent  screen  was  within  8%.  Using  the  developed  simulation 
technology,  we  will  be  able  to  produce  simulated  mammograms  or  other  x-ray  images 
that  can  be  used  to  evaluate  CAD  methods  and  techniques. 

Keywords:  x-ray  imaging,  mammographic  imaging,  fluorescent  screens,  computer  aided 
diagnosis,  computer  simulation 

1.  Introduction: 

Many  research  groups  have  been  investigating  CAD  methods  to  help  radiologists 
improve  their  detection  and  diagnosis  of  breast  cancer  from  mammograms  and  various 
methods  and  technologies  have  been  proposed  [1—5].  A  limitation  to  these  quantitative 
mammogram  analysis  techniques  is  a  lack  of  the  set  of  standardized  mammograms  to 
evaluate  the  methods  developed  from  different  research  labs.  M.  L.  Giger  [4]  reviewed 
the  computer-aided  diagnosis  in  breast  imaging  and  emphasized  the  importance  to 
examine  the  performance  of  the  CAD  system  and  evaluate  their  ultimate  contribution  to 
the  breast  cancer  diagnosis  before  using  them  in  clinical  practice.  Although  developing  a 
common  database  of  mammograms  is  helpful  in  solving  the  problem,  the  data  set  will 
still  be  less  than  ideal  because  the  lack  of  accurate  truth  information  (exact  size,  shape, 
and  location)  for  the  lesions  in  the  mammograms.  So,  computer  simulation  technology  is 
considered.  In  this  paper,  we  propose  a  technology  to  numerically  simulate  the  physics 
process  of  a  mammographic  imaging  system.  This  technology  can  be  used  to  generate 
simulated  mammograms  with  truth  information  of  lesions  to  evaluate  CAD  methods 
developed  in  labs. 

The  physics  process  we  simulated  is  the  conversion  process  from  x-ray  image  input 
to  visible  light  image  output  in  most  mammographic  imaging  systems.  In  clinical  practice, 
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most  mammograms  are  taken  from  a  system  that  employs  a  fluorescent  intensifying 
screen,  as  shown  by  Figure  1.  In  these  systems,  the  x-ray  quanta  transmitted  through  the 
object,  for  instance,  breast,  impinge  upon  the  fluorescent  intensifying  screen.  The  x-ray 
quanta  that  interact  in  the  screen  will  excite  the  electrons  from  the  valence  band  to  the 
conduction  band  inside  the  screen.  In  an  efficient  screen,  many  of  those  electrons  will 
emit  visible  light  when  returning  to  valence  band.  The  fluorescent  intensifying  screen 
works  as  a  linear  transducer  that  converts  absorbed  x-ray  quanta  into  light  photons.  The 
intensifying  screen  shown  by  Figure  1  reduces  the  patient  radiation  dose  compared  to 
using  x-ray  film  alone.  However,  it  also  has  negative  affect.  First,  the  signal  to  noise  ratio 
(SNR)  of  final  x-ray  images  may  be  reduced  slightly  because  the  x-ray  quantum  noise 
may  be  transferred  more  efficiently  than  signals  [6-8],  Secondly,  the  screen  can  also 
bring  some  special  noise  such  as  the  screen-structure  noise,  secondary-quantum  noise,  etc 
[9,  10].  Thirdly,  the  scatter  of  the  light  photons  created  in  the  screen  reduces  spatial 
resolution  of  the  final  images.  All  these  influences  from  the  fluorescent  intensifying 
screen  are  critical  to  diagnose  small  lesions  at  an  early  stage.  So,  it  is  important  to 
simulate  the  physics  processes  inside  the  fluorescent  intensifying  screen  accurately  if 
using  computer  simulated  images  as  the  data  to  evaluate  CAD  technologies. 

In  this  paper,  we  describe  a  method  to  numerically  simulate  the  physics  processes 
inside  the  fluorescent  intensifying  screen  used  in  mammography.  The  method  can  also  be 
used  to  simulate  other  x-ray  imaging  system  using  the  x  rays  less  than  40  keV.  We 
developed  software  using  this  technology  to  simulate  the  conversion  of  x-ray  quantum  to 
visible  light  photons  inside  the  fluorescent  intensifying  screen.  This  work  was  based  on 
the  model  of  the  spatial-frequency-dependent  detective  quantum  efficiency  of  phosphor 
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screens  that  we  had  previously  [9,  10].  The  technology  can  be  used  to  generate  the 
simulated  x-ray  image  with  accurate  truth  information  for  lesions. 

2.  Simulation  Model 

The  model  of  the  spatial-frequency-dependent  detective  quantum  efficiency  of 
phosphor  screens  uses  six  stages  to  describe  the  physics  processes  inside  a  fluorescent 
intensifying  screen  [3].  The  physics  processes  of  converting  x-ray  quanta  to  light  photons 
inside  a  fluorescent  screen  are  separated  into  two  elementary  processes:  amplification 
process  and  scattering  process  [11].  Stages  1  to  4  form  a  quantum  number  amplification 
process  and  stages  5  and  6  form  a  photon  scattering  process.  Because  we  limited  our 
simulation  to  mammographic  imaging  system  or  x-ray  imaging  systems  using  the  x  rays 
less  than  40  keV,  the  K  or  L  fluorescence  of  the  screens  can  be  ignored. 

The  first  stage  is  where  x-ray  quanta  are  incident  on  a  fluorescent  intensifying  screen 
of  thickness  T.  We  use  u(  ri  . ,  E)  to  express  the  distribution  of  x-ray  fluence  in  this  stage, 

where  rt  .  is  a  two-dimensional  vector  describing  a  pixel  position  when  x-ray  quanta  with 

energy  E  are  incident  on  the  screen.  If  the  pixel  size  is  A  a  in  both  width  and  height  and 
using  i  as  the  horizontal  index  and  j  as  the  vertical  index,  we  have 

r.  .  =  A  ax  ii  +  A  ax  jj  .  i  is  the  unit  vector  in  horizontal  direction  and  j  is  the  unit  vector 
in  vertical  direction.  The  second  stage  is  the  fraction  of  the  input  quanta  that  interact  in 
the  fluorescent  screen.  We  use  v(f; } ,  E)  to  express  the  distribution  of  the  x-ray  quanta  that 

will  interact  in  the  screen  after  entering  the  screen  at  a  pixel  area  Aa .  The  third  stage  is 
where  x-ray  quanta  will  interact  at  different  depths,  t,  inside  the  fluorescent  screen.  We 
use  }  ,t,  E)  to  express  the  distribution  of  x-ray  quanta  in  a  voxel  AaAt ,  centered  at 
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(rUit) ,  at  this  stage.  The  fourth  stage  is  where  the  absorbed  x-ray  energy  is  partially 
converted  into  light  photons.  We  use  ;.,t)to  express  the  distribution  of  the  light 

photons  created  at  this  stage.  From  stage  1  to  stage  4,  the  input  x-ray  quanta  become 
visible  light  photons.  Because  the  energy  of  light  photons  is  much  smaller  than  that  of  x- 
ray  quanta,  the  number  of  visible  light  photons  is  much  larger  than  the  number  of  primary 
input  x-ray  quanta.  So,  stage  1  to  stage  4  forms  a  quantum  number  amplification  process. 
u{rtj ,  E )  is  the  input  of  this  quantum  number  amplification  process  and  xty  ,,t)  is  the 

output  of  this  quantum  number  amplification  process. 

The  fifth  stage  is  the  fraction  of  light  photons  created  at  each  voxel  that  escape  to  the 
output  surface  of  the  intensifying  screen.  The  photon  distribution  after  this  escape  can  be 
expressed  as  y(jrr  j,t) .  The  sixth  stage  is  where  the  photons  escaping  spread  over  the 
output  and  give  the  final  light  output.  The  distribution  of  photons  can  be  expressed 
as  zir-j ) ,  where  ?:  J  is  the  position  vector  of  an  output  light  image  pixel.  These  two  stages 

form  the  scattering  process  of  the  light  photons  created  inside  the  screen  and  characterize 
scattering  and  absorption  of  the  light  photons  created  from  the  energy  of  input  x-ray 
quanta  as  shown  by  Figure  2.  Swank  [12]  derived  a  formula  to  calculate  an  1- 
dimensional  modulation  transfer  function  [MTF(/)]  for  modeling  the  diffusion  and 
absorption  process  of  light  photons  inside  phosphor  by  assuming  an  1 -dimensional  (line) 
input  source.  To  simulate  the  2-dimensional  mammograms  or  other  x-ray  images,  a 
formula  to  modeling  2-dimensional  modulation  transfer  function  of  diffusion  and 
absorption  of  light  photon  is  required.  We  derived  this  formula,  Equation  (1),  using  the 
same  method  as  Swank  used  (please  see  Appendix  A).  In  fact,  Equation  (1)  is  the 
solution  of  Boltzmann  differential  equation  in  spatial  frequency  domain. 


6 


(puv{u,v,t)  = 


zp\  [(0)  +  zp0  )e°*  +  (CO -  zp0  )e~ox  ] 

{to + Tp0 \o) + zpx  )e6Sr  -  {(d-zps )(C0 -  zpx )e~ 


where  (pu  v  ( u,v,t )  is  the  Fourier  Transformation  of  the  product  of  light  photon  density 


and  diffusion  coefficient  (see  Appendix  A;)  u  and  v  are  spatial  frequency;  t  is  the  depth 
inside  a  fluorescent  screen;  z  is  the  inverse  relaxation  length  that  describes  the  probability 
of  optical  absorption  occurring;  ry  is  related  to  spatial  frequency  and 
O)1  =  a2  +4 7Z(u2  +  v2)  where  <7  is  the  reciprocal  diffusion  length  and  characterizes 


optical  scattering;  T  is  the  thickness  of  the  screen;  p{= - -  where  r;  is  the  reflectivity  of 

1  +  r, 

the  boundary  of  the  screen;  ( i  =  0 )  represents  the  interface  closest  to  the  x-ray  source  of 
the  screen  and  ( i  =  1)  represents  the  output  interface. 

Equation  (1)  can  be  used  to  simulate  scattering  and  absorption  process  of  light 
photons  created  inside  a  fluorescent  intensifying  screen.  Parameters  z,  p0 ,  and  a 
depend  on  the  type  of  the  phosphor,  the  size  and  the  shape  of  the  phosphor  crystals,  the 
absence  or  addition  of  some  light  absorbing  materials,  and  the  reflectivity  of  input  and 
output  interfaces  of  the  screen.  Those  parameters  characterize  the  screens  used  in 
mammographic  imaging  systems  or  other  x-ray  imaging  systems. 

Figure  3  shows  two  examples  of (pu  v(u,v,t) function  at  different  screen  depths.  The 


parameters  to  calculate  Figure  3  are  listed  in  Table  1.  They  are  chosen  to  simulate  a 
Mini  R  2000  intensifying  screen.  The  left  panel  of  Figure  3  was  calculated  at  screen 
depth  t  =  \0pm  and  the  right  panel  of  Figure  3  was  calculated  at  screen  depth  t  =  80 pm  . 
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Table  1  Simulation  Parameters 


X-ray  Energy  E 

20.0  KeV 

Pixel  Size 

17x17  pm 2 

Image  Size 

256x256  pixel 

Linear  Attenuation  p 

37.7  1/mm 

Conversion  Coefficient  £ 

0.15 

Screen  Thickness  T 

85  pm 

Reciprocal  Diffusion  Length  a 

11  1/mm 

Inverse  Relaxation  Length  t 

50  1/mm 

Input  Ratio  of  Reflectivity  pQ 

0.90 

Output  Ratio  of  Reflectivity  px 

0.70 

Sub-layer  Thickness 

1  pm 

3.  Simulation  Algorithm 

Our  algorithm  to  simulate  the  physics  processes  inside  a  fluorescent  intensifying 
screen  is  to  conceptually  divide  the  fluorescent  screen  into  multiple  thin  sub-layers  shown 
by  Figure  4.  Each  sub-layer  is  thin  enough  to  be  of  a  fixed  depth  t.  The  closer  to  the 
input  surface  of  the  screen  the  sub-layer  is,  the  higher  the  input  x-ray  fluence  to  the  sub¬ 
layer  is.  The  diffusion  and  absorption  of  light  photons  created  at  each  sub-layer  depend 
on  how  far  the  sub-layer  is  from  the  output  surface  of  the  screen.  There  will  be  more  light 
photons  diffused  and  absorbed  if  the  sub-layer  is  far  from  the  output  surface  of  the  screen. 
We  used  five  cascading  mathematical  transformations  to  simulate  the  physics  processes 
of  converting  x-ray  quantum  input  to  the  final  light  image  output  at  each  sub-layer.  The 
output  light  fluence  of  the  whole  intensifying  screen  was  considered  as  the  accumulation 
of  the  light  output  of  every  sub-layer.  The  simulation  software  integrates  the  light  outputs 
of  all  sub-layers  into  a  final  output  image  of  the  phosphor  screen. 
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The  five  cascading  transformations  used  are  x-ray  interaction  mathematical 
transformation,  x-ray  quantum  to  light  photon  transformation,  Fourier  transformation, 
light  photon  escape  and  scattering  transformation,  and  inverse  Fourier  transformation. 

The  first  transformation  used  was  x-ray  interaction  mathematical  transformation.  It 
simulates  the  physics  process  from  stage  1  to  stage  3  of  the  model  of  the  spatial- 
frequency-dependent  detective  quantum  efficiency  of  phosphor  screens  [9,  10].  If  the 
number  of  x-ray  quanta  is  N0  before  the  x-ray  beam  is  incident  to  the  voxel  centered  at 

(. x ,  y)  of  a  sub-layer  thick  At  at  depth  t,  according  to  linear  attenuation  theory  [13],  the 
number  of  x-ray  quanta  penetrating  the  voxel  is 
N  =  iV0e“M'  (2) 

where  N  is  the  number  of  the  x-ray  quanta  penetrating  from  this  voxel;  fl  is  the  linear 
attenuation  coefficient  of  the  intensifying  screen;  At  is  the  thickness  of  the  sub-layer 
concerned.  From  Equation  (2),  we  can  get  the  number  of  the  x-ray  quanta  removed  from 
the  beam  in  the  voxel  observed: 

AN  =  juNAt  =  /uN^Ate'^  (3) 


where  AN  is  the  number  of  the  x-ray  quanta  removed  from  the  beam  after  it  incidents  to 
the  voxel  centered  at  {x,  y )  of  the  sub-layer  at  depth  t  of  the  screen. 

From  Equation  (3),  we  can  get  the  fraction,  PremovingX-r„y ,  that  the  x-ray  quanta  are 

removed  from  the  beam  at  the  voxel  centered  at  (x,  y)  and  at  depth  t  of  the  intensifying 
screen.  It  can  be  expressed  as 


P  .  x  =  —  =  juAte-^ 

removing*  -ray  j.r  r”—* 

^ 0 


(4) 
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Using  this  PremovingX-ray » the  x-ray  interaction  mathematical  transformation  can  be 
expressed  as  Equation  (5). 

w(x,  y,  E,  t )  =  fjAte^' u(x,  y,  E )  (5) 

The  second  mathematical  transformation  was  x-ray  quantum  to  light  photon 
transformation.  It  simulates  the  creation  of  the  light  photons  from  the  absorbed  energy  of 
x-ray  quanta.  According  to  the  assumptions  of  the  model  of  spatial-frequency-dependent 
detective  quantum  efficiency  of  phosphor  screens,  every  photon  created  inside  the  screen 
has  energy  2.4  eV  or  2.4xl0~3  keV  [10,  14].  Then,  the  number  of  photons  created  by 
absorbing  one  x-ray  quantum  with  energy  E  is 

E 

q(x,  y,t)=  £w(x,  y,  E,  t )  (6) 

2.4x10 

where  q(x,  y,  t)  is  the  distribution  of  photons  created  by  absorbing  x-ray  quanta;  E  is  the 
energy  of  the  x-ray  quantum  in  keV;  £  is  the  conversion  coefficient  of  the  screen  [10]. 

Mathematical  transformations  3  to  5  simulate  the  diffusion  and  absorption  of  the 
photons  created  inside  the  intensifying  screen.  These  photons  need  to  escape  from 
absorption,  diffusion,  and  reflection  before  contributing  to  the  output  light  image.  We  can 
use  Equation  (1)  to  simulate  the  affects  of  this  escape.  However,  Equation  (1)  is  the 
solution  of  the  Boltzmann  diffusion  equation  in  spatial  frequency  space  (please  see 
Appendix  A).  We  either  use  Equation  (1)  as  a  kernel  to  do  convolution  with  the  output  of 
transformation  2  or  transform  the  output  of  transformation  2  into  spatial  frequency  space. 
To  save  the  computation  time  from  convolution,  we  transform  the  output  of  the  second 
transformation  into  spatial  frequency  space.  We  used  standard  2-demisonal  FFT 
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technology  to  do  this  transformation  [15].  So,  our  third  mathematical  transformation  used 
was  a  Fourier  transformation  shown  by  Equation  (7) 

Q(u,v,t)  =  F{q(x,y,t)}  (7) 

where  u  is  the  spatial  frequency  in  horizontal  direction;  v  is  the  spatial  frequency  in 
vertical  direction,  t  is  the  depth  position  of  the  sub-layer. 

The  fourth  mathematical  transformation  used  was  in  spatial  frequency  domain  and 
written  as: 

Z(u,v,t)  =  Q(u,v,t)<puv(u,v,t)  (8) 

where  Z(u,  v,  t)  is  the  Fourier  Transformation  of  the  final  output  image;  <pu  v(u,v,t )  is 
defined  in  Equation  (1). 

The  fifth  mathematical  transformation  used  was  an  Inverse  Fourier  Transformation 
shown  by  Equation  (9). 

z(x',y',t)  =  F~l{Z(u,v,t)}  (9) 

where  z(x',  y',t )  is  the  final  output  image  from  the  concerned  sub-layer  at  depth  t  of  the 
intensifying  screen.;  (x\  y')  is  the  pixel  position  on  the  output  interface  of  the 
intensifying  screen.  Figure  5  shows  the  five  cascading  mathematical  transformations  in 
our  simulation  software. 

Figure  6  shows  the  data  flow  of  our  x-ray  imaging  simulation  software.  Users  can 
set  characteristics  of  the  fluorescent  intensifying  screen,  select  x-ray  inputs,  and  assign 
the  thickness  of  the  sub-layer  through  user  interface.  After  loading  the  x-ray  input  data 
and  initializing  the  simulation  system  according  to  user’s  selection,  the  simulation 
software  transfers  the  pixel  values  of  the  input  x-ray  image  into  the  number  of  x-ray 
quanta.  Then,  the  simulation  kernel  comes  into  the  computation  loops  of  sub-layers.  In 
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each  loop,  the  simulation  software  first  sets  the  thickness  of  the  sub-layer  and  its  depth 
position  in  the  intensifying  screen.  Then,  the  simulation  software  performs  the  five 
mathematical  transformations  to  the  input  x-ray  image  data.  Finally,  the  simulation 
software  accumulates  the  output  of  each  sub-layer  to  the  output  buffer.  After  processing 
the  current  sub-layer,  the  simulation  software  checks  if  there  are  more  sub-layers  required 
to  be  dealt  with.  If  there  are  untreated  sub-layers,  the  simulation  software  comes  into  the 
next  loop  of  sub-layer  computation.  Otherwise,  the  program  goes  into  the  output 
operation.  The  output  operation  includes  transforming  the  number  of  light  photons  into 
optical  densities,  buffering  the  data,  etc. 

4.  Simulation  Input  and  Output 

The  x-ray  input  of  the  simulation  software  is  obtained  by  calibrating  the  directly 
contacting  film  x-ray  images.  We  first  made  a  step  image  using  directly  contacting  x-ray 
film  to  calibrate  the  x-ray  exposure  and  pixel  values.  The  exposure  for  every  step  was 
recorded  when  we  made  the  step  image.  After  developing  the  step  image  film,  we 
measured  the  optical  density  of  the  step  image  film  and  digitized  the  step  image  film  into 
a  12  bit  4000x5000  digital  image.  On  this  digital  image,  every  step  corresponds  to  a  pixel 
value.  By  the  exposure  records,  we  built  the  relationship  between  pixel  values  and  x-ray 
exposures  as  shown  by  Table  2.  The  smallest  pixel  value  is  0  and  the  biggest  pixel  zero  is 
4095.  According  to  the  relationship  between  x-ray  exposure  and  x-ray  quantum  number, 
we  get  the  x-ray  quantum  number  for  every  special  pixel  value  of  our  digitalized  step 
image  [16],  From  reference  [17]  and  our  own  test,  we  used  4.4x10"  -6.9x10" 
quanta  /(mm2 R)  at  20  keV  to  transform  the  x-ray  exposure  into  x-ray  quantum  number. 
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Table  2  X-ray  exposure  to  pixel  value 


Exposure  (  R  ) 

Pixel  Value 

Step  1 

1.93 

3570 

Step  2 

1.07 

3370 

Step  3 

0.68 

2984 

Step  4 

0.47 

2216 

During  buffering  the  input  x-ray  image  in  simulation,  the  simulation  software  translates 
the  pixel  values  of  the  selected  x-ray  digital  image  based  on  the  calibration,  using 
interpolating  and  extrapolating  technology.  When  doing  interpolation  and  extrapolation, 
we  assume  the  x-ray  exposure  information  is  in  the  domain  of  the  linear  portion  of  the 
film’s  characteristic  curve  [17]. 

The  simulation  output  image  was  obtained  by  transforming  the  number  of  light 
photons  of  each  pixel  into  optical  density  using  curve  of  optical  density  to  photon  fluence 
shown  by  Figure  7,  taken  from  reference  [17].  Then,  the  optical  density  was  transformed 
into  pixel  value  using  the  curve  of  the  pixel  value  to  optical  density  shown  by  Figure  8. 
This  curve  was  obtained  by  calibrating  a  step  x-ray  image  from  an  imaging  system 
employing  an  intensifying  screen. 

5.  Simulation  Assurance 

We  used  three  methods  to  check  our  simulation  algorithm  and  the  simulation  software. 
The  three  testing  methods  are:  using  known  images,  checking  the  modulation  transfer 
function  of  the  simulation  software  [18-21],  and  checking  the  noise  power  spectrum  [22- 
26]  of  the  simulation  software.  To  compare  with  the  measured  modulation  transfer 
function  and  noise  power  spectrum  from  Mini  R  2000  screen-film  system,  we  set  the 
simulation  parameters  as  listed  in  Table  1. 

Figure  9  and  Figure  10  show  the  comparison  of  the  simulated  image  with  the  x-ray 
image  taken  from  a  Mini  R  2000  screen-film  mammographic  imaging  system.  The  left 
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panel  of  Figure  9  is  a  phantom  image  taken  from  the  real  screen-film  system.  The  right 
panel  is  a  simulated  image  of  the  same  phantom  object.  The  two  images  look  similar.  The 
left  panel  of  Figure  10  is  the  display  values  from  a  row  of  a  screen-film  phantom  image. 
The  width  of  its  pixel  is  50  /An .  The  right  panel  of  Figure  10  is  the  display  values  from  a 
row  of  a  simulated  image  of  the  same  phantom  object.  Because  the  width  of  the  simulated 
image  is  17  /An ,  the  display  values  plotted  are  average  per  3  pixels.  The  comparison 
shows  that  the  simulated  image  has  same  imaging  futures  as  the  image  from  the  real 
screen-film  mammographic  imaging  system. 

To  measure  the  modulation  transfer  function  (MTF)  of  our  simulation  operation,  we 
first  generated  a  point  image  using  the  computer.  This  image  has  only  one  bright  pixel  at 
its  center  with  the  maximum  pixel  value  4095  which  is  corresponding  to  x-ray  exposure 
4. 1 8R  according  to  Table  2.  This  exposure  is  equal  to  about  4600  x-ray  quanta  if  the 
pixel  size  is  17/Anxl7/An  and  the  x-ray  fluence  is  6.9x10"  quanta  /(mm2 R)  at  20  keV. 
We  ran  the  simulation  software  using  this  point  image  as  the  x-ray  input.  Then,  we  did  a 
Fourier  Transformation  to  the  simulation  output  of  this  point  image.  According  to  the 
definition  of  modulation  transfer  function,  this  Fourier  Transformation  is  the  modulation 
transfer  function  of  our  simulation  software.  Figure  11  shows  the  input  point  image  and 
the  Fourier  Transformation  of  the  simulation  output.  The  left  panel  of  Figure  11  is  the 
point  image  and  the  right  panel  is  the  Fourier  transformation  of  the  simulation  output  of 
the  point  image.  Figure  12  shows  the  modulation  transfer  factor  of  our  simulation 
software,  initialized  with  the  parameters  listed  in  Table  1,  and  the  modulation  transfer 
factor  measured  from  a  Mini  R  2000  intensifying  screen.  By  comparing  the  modulation 
transfer  function  of  simulation  software  with  the  measured  data,  we  can  see  that  the 
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simulation  software  works  correctly.  The  difference  between  the  modulation  transfer 
factor  of  our  simulation  software  and  the  measured  values  is  within  4%. 

To  obtain  the  quantum  noise  power  spectrum  of  the  simulation  software,  we  created  a 
noise  image  with  the  normal  random  number  generator  of  MATLAB  as  the  input  of  our 
simulation  operation.  We  subtracted  the  mean  from  the  simulation  output  and  squared  the 
result.  Then,  we  did  Fourier  Transformation  and  divided  the  result  by  the  maximum  value 
in  spatial  frequency  domain.  Figure  13  shows  one  of  our  noise  input  images  and  the 
simulation  outputs.  The  left  panel  of  Figure  13  is  the  input  noise  image  and  the  right 
panel  is  the  output  of  the  simulation.  We  can  see  the  blurring  resulting  from  the  simulated 
conversion  process  from  x-ray  quanta  to  light  photons  and  the  scattering  of  light  photons 
inside  the  simulated  intensifying  screen. 

Figure  14  shows  the  relative  noise  power  spectrum  of  our  simulation  software.  The 
left  panel  is  an  original  noise  power  spectrum  of  our  simulation  software  from  the 
measurement.  The  right  panel  of  Figure  14  shows  the  smoothed  relative  noise  power 
spectrum  of  the  simulation  software  and  the  measured  relative  noise  power  spectrum  of  a 
real  screen.  The  measured  data  are  from  Kodak  Min-R  2000  screen  [23].  The  smoothed 
relative  noise  power  spectrum  of  simulation  software  is  calculated  by  two  steps  of 
averaging  computation.  The  first  step  is  that  we  put  4  consecutive  data  of  a  relative  noise 
spectrum  into  one  group  and  calculated  the  average  of  every  group.  The  second  step  is 
that  we  used  1 1  random  noise  input  images  as  shown  by  the  left  panel  of  Figure  13  and 
obtained  1 1  relative  noise  power  spectra  of  simulation  software.  Then,  we  averaged  the 
1 1  spectra.  We  can  see  that  the  difference  between  the  relative  noise  power  spectrum  of 


our  simulation  software  and  the  measured  values  is  within  8%. 
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6.  Summary 

This  paper  presented  a  method  to  simulate  the  physics  processes  of  converting  x-ray 
quantum  image  to  visible  light  photon  image  inside  a  fluorescent  intensifying  screen 
based  on  previous  work  in  modeling  the  quantum  efficiency  of  the  fluorescent 
intensifying  screen.  We  divided  the  fluorescent  intensifying  screen  into  multiple  thinner 
sub-layers  and  used  five  mathematical  transformations  to  simulate  the  physics  process  of 
every  sub-layer  that  converts  the  x-ray  input  quanta  to  light  photons  and  the  process  of 
the  diffusion,  absorption,  and  reflectivity  of  the  light  photons.  We  integrated  the  light 
photons  of  every  sub-layer  over  the  whole  screen  into  the  final  output  image.  Using  the 
algorithm  of  this  paper,  we  developed  the  simulation  software  in  Linux  operating  system, 
using  C++.  Three  methods  were  used  to  check  our  simulation  technology.  Qualitatively, 
the  simulated  and  actual  images  of  test  objects  looked  similar.  Quantitatively,  the 
difference  between  our  simulation  software  and  the  measured  values  of  a  real  intensifying 
screen  was  within  4%  in  modulation  transfer  function  and  within  8%  in  relative  noise 
power  spectrum. 

Based  on  this  x-ray  imaging  system  simulation  technology,  it  is  possible  to  generate 
some  simulated  mammograms  or  other  x-ray  images  with  accurate  truth  information  for 
lesions.  The  simulated  x-ray  images  can  be  used  to  evaluate  computer  aided  diagnose 
(CAD)  technology.  The  further  work  includes  modeling  the  intensifying  screen  more 
accurately,  measuring  the  linear  attenuation  coefficient  of  the  malignant  lesions  and 
benign  lesions  accurately,  and  researching  the  algorithms  to  add  the  lesions  to  x-ray 
images  for  x-ray  imaging  simulation  software.  We  also  need  to  develop  some  algorithm 
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to  simulate  the  more  serious  random  noise  of  x-ray  quanta  in  screen-film  system  than  that 
of  directly  contacting  film  images. 

7.  Acknowledgement 

This  work  was  funded  in  part  by  grant  from  the  US  Army  DAMD 17-99-9 122  and  the 
NIH  T32  CA09649. 

Appendix  A:  Derivation  of  light  photon  diffusion  and  absorption  model 

The  Boltzmann  differential  equation  of  photons  created  inside  a  screen  can  be  written 
as  Equation  (Al)  if  the  observed  points  are  relatively  far  from  the  source  ( i.e .  the  distance 
between  the  observed  point  and  source  is  larger  than  a  mean  free  path)  and  the  observed 
point  is  not  close  to  the  boundary  [27]. 

V2A0(F)--^A0+lso=0  (Al) 

where  A0is  the  particle  density;  S0is  an  isotropic  source  of  particles;  r  is  the  position 

vector  of  the  observed  point;  L  is  the  diffusion  length;  D  is  the  diffusion  coefficient. 

Both  L  and  D  are  related  to  the  transport  mean  free  path  of  the  particles. 

Inside  a  fluorescent  screen  (see  Figure  2),  the  Boltzmann  equation  can  be  simplified 
into  Equation  (A2)  by  combining  L  and  D  into  one  coefficient:  reciprocal  diffusion  length. 
-V2  </>(?)  +  a2  0(r)  =  S(r)  (A2) 

where  0(r)  is  the  product  that  the  light  photon  density  times  the  diffusion  coefficient; 


S(r)  is  the  source  function;  a  is  the  reciprocal  diffusion  length  and  o  ■ 


D 

L 1 


If  the  input  image  is  u(x,  y) ,  the  sampling  of  the  input  image  can  be  expressed  as 
u  +  (x,  y)  =  u(x,  y)^  8(x-ia,y-  ja )  =  ^  uL .  (ia,  ja)8(x  -  ia,  y  -  ia)  (A3) 

‘•j  Uj 
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where  u  +  (x,  y )  is  a  sequence  of  scaled  8  functions.  So,  we  only  need  to  derive  a 
solution  of  Equation  (A2)  when  the  source  is  a  2-dimensional  8  function. 

If  the  source  in  Equation  (A2)  is  a  point  source  andS(r)  =  8(x)S(y)8(t  -t0).  Equation 
(A2)  becomes 

-  V  V  +  c2( i>  =  8(x,  y,t-tQ)  =  8(x)8(y)8(t  -t0)  = 

T  1  r  ■  Y  1  r  •  (A4) 

In  Equation  (A4),  we  use  a  definition  of  8  function,  Equation  (A5),  to  substitute  8(x) 
and  8(y)  [28], 

=  -L  \ei2™du  (A5) 

2k  j 

If  we  assume 

0(x,  y,t)  =  ~T  |K,  (t)ei2*ULX+v>) dudv ,  (A6) 

and  substitute  0(x,  y,t)  of  Equation  (A4)  with  Equation  (A6),  we  can  get  a  new  equation, 
Equation  (A7). (Please  see  Appendix  B) 

-  V^„,v (0  +  [4^2 (u2+v2)  +  a2 }(pu  v ( t )  =  8(t  - 10 )  (A7) 

Equation  (A7)  is  Equation  (A4)  in  spatial  frequency  space.  Using  the  same  method  as 
Swank  used  [12],  we  can  obtain  Equation  (1)  as  the  solution  of  Equation  (A7). 


Appendix  B:  Derivation  of  Equation  (A7) 

Using  Equation  (A6),  we  can  get  Equation  (Bl)  for  operator  V 2 . 
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V^Oc, y,t)  =  V2]-^ 

;  (B1) 

In  Equation  (Bl),  we  exchanged  the  order  of  derivative  and  integration. 

Using  Equation  (Bl)  and  Equation  (A5),  Equation  (A4)  can  be  rewritten  into 


- ^7  fjv2k.v(0e— } \ludv  +  rx2~  \\<pu,M)ei2My)dudv 

L_  ^e^^dudv\s(jt-to) 


(B2) 


9  _  ^  _  0  _ 

Considering  V  =  —i  H - j  H — A: ,  the  derivative  in  Equation  (Bl)  can  be  calculated 

ox  ay  at 


by  Equation  (B3) 

vfeM,(  t)ei2*^) 

=  (puv(t){i2m)ei2*{m+vy)T  +  ^,v(0(/2^v)e'' 


1 27T(/c.r+iry)  7  i  3%.v(0 


;  +■ 


dt 


2/r(i«+v>’)  £ 


(B3) 


Using  Equation  (B3),  we  can  get 

v2fc,v(Oe,2*(“+^) 

=  V  j  <puv  (t){i2m  )ei2*(ux+vy)i  +  <puv  (t)(i2nv)eiW"x+vy)  j  +  d(p^(*t)  ei2*^y>  i 


i2ff(ux+vy) 


=  <P»,V  (0[-  4^t2m2  \n*^  +  (puv  (r)[-  4n2v2\' 

=  -4# 2 (u2  +  v2  }pil  v (t)e,2aUlx+vy)  +  (y2<pu  v (t)]en*{ux+vy) 


+ 


dt 

d2<pu,.(t) 
dt2 


27T(lLX+vy) 


(B4) 


7  d2 

In  Equation  (B4),  we  defined  Vr  =  — - . 

dt 


Using  Equation  (Bl)  and  Equation  (B4),  we  can  rewrite  Equation  (A4)  into: 
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- l-jj{-4x2(u2  +v2kv(0ei2"(“+'7)  +(v>HV(oV2"M}/Hrfv 

4;r  .  r  i  1  <B5) 

+  rr2  j J^„,v(0«i2,r(“+’7)^v  =  |  jjeimux+vy)dudv  \8{t  - 10 ) 

Combining  the  integration  on  the  left  of  Equation  (B5)  and  writing  8{t  - 10 )  into  the 
integration  on  the  right  of  Equation  (B5),  Equation  (B5)  becomes: 

~Y  jf{-  V>„  v (0  +  [4^2 (u2  +v2)  +  a2 }<pu  v 

4;r  (B6) 

=  -k  f  fa*  - 10  )eim,,x+vy)dudv 
An 

From  Equation  (B6),  we  can  get  Equation  (A7)  or  (B7)  with  (O2  =  An2(u 2  +  v2)  +  o2 

-  v,V„,v  (0 + ®V,,v  (0  =  #(t-t0)  (B7> 
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Phosphor  Screen 


Figure  1:  The  x-ray  imaging  system  employing  a  fluorescent  intensifying  screen. 


Figure  2  Diffusion  of  photons  created  by  x-ray  quanta.  The  photons  created  at 
O,  y,  t )  may  diffuse  to  (V,  y ',0) . 
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p(uv)  et  1 Q  i»n  h  a  es>»n  screen  f(u.v)  ai  80>r>  In  a  65  yn  screen 


Figure  3  Photon  diffusion  and  absorption  function  inside  a  fluorescent  intensifying 
screen  with  Equation  (1).  The  parameters  of  the  screen  are  listed  in  Table  1.  The  left 
panel  was  calculated  when  t  =  10 /im  and  the  right  panel  was  calculated  when?  =  80 flm . 
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Figure  4  The  fluorescent  intensifying  screen  is  conceptually  divided  into  multiple  thin 
sub-layers.  All  the  sub-layers  are  independent  of  each  other.  Every  sub-layer  has  different 
input  x-ray  fluence,  absorbs  different  number  of  x-ray  quanta,  and  outputs  different  light 
fluence.  The  diffusion  and  absorption  of  light  photons  of  a  sub-layer  depend  on  how  far 
the  sub-layer  from  the  output  interface  of  the  screen. 


Input  Image 


Figure  5  Five  mathematical  transformations  used  to  simulate  the  processes  of  converting 
input  x-ray  quanta  to  light  photons.  The  input  x-ray  image  to  every  sub-layer  will  pass 
these  five  cascading  mathematical  transformations  to  contribute  to  the  final  output  light 


image. 
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Figure  6  The  data  flow  of  the  simulation.  After  buffering  the  x-ray  input  image  and 
initializing  the  simulation  system  according  to  user’s  assignments,  the  simulation 
software  goes  into  the  loop  to  do  the  computation  for  every  sub-layer  of  the  intensifying 
screen.  The  final  light  image  output  of  the  whole  intensifying  screen  is  the  accumulated 
output  of  all  sub-layers. 
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The  Phantom  Image  from  A  Mini  R2000  System  The  Simulated  Image  of  The  Same  Object 


Figure  9  Comparison  of  the  simulated  phantom  image  with  the  image  taken  from  a 
screen-film  system.  The  left  panel  is  the  image  taken  from  a  Mini  R  2000  screen-film 
system.  The  right  panel  is  a  simulated  image. 


33 


Screen-film  Phantom  Image  Display  Value  Simulated  Phantom  Image  Display  Value 


Figure  10  Comparison  of  the  display  values  between  the  screen-film  phantom  image  and 
simulated  phantom  image.  The  left  panel  is  the  display  values  from  a  row  of  screen-film 
phantom  image.  The  right  panel  is  the  display  value  from  a  row  of  simulated  image  of  the 
same  phantom  object. 
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Computer  generated  poH  image  The  Fourier  Transformation  of  the  skmttlon  oiiptl 


Figure  11  the  measurement  of  the  modulation  transfer  function  of  the  simulation 
software.  The  left  panel  is  the  input  image  used  to  measure  the  modulation  transfer 
function.  The  right  panel  is  the  Fourier  transformation  of  the  simulation  output,  the 
modulation  transfer  function  of  the  simulation  software. 
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Figure  12  Comparison  of  the  modulation  transfer  factor  of  the  simulation  software  with 
the  measured  modulation  transfer  factor.  The  solid  curve  is  the  modulation  transfer  factor 
of  the  simulation  software,  [x]  is  the  measured  modulation  transfer  factor  from  a  Min-R 
2000  screen[23]. 
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A  2D-Noise  Spectrum  of  The  Simulation  System  Smoothed  Relative  Noise  Spectrum  of  The  Simulation  System 


Figure  14  The  noise  power  spectrum.  The  left  panel  is  a  2D  noise  power  spectrum  of  the 
simulation  software.  The  right  panel  compares  the  smoothed  relative  noise  power 
spectrum  of  the  simulation  software  with  the  measured  relative  noise  power  spectrum 
from  a  Mini  R2000  intensifying  screen.  [  A  ]  is  the  relative  noise  power  spectrum  of  the 
simulation  software,  [x]  is  the  relative  noise  spectrum  measured  from  a  Min-R  2000 
screen  [23]. 
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3D  Modeling  and  Simulation  of  Characteristic  x-ray  re-absorption  by  Monte 

Carlo  Simulation 

Yinghui  Zhang  PhD.  and  Robert  M.  Nishikawa  PhD 
The  radiology  department,  the  University  of  Chicago,  IL  60637 

Abstract 

Many  x-ray  detectors  use  a  phosphor  screen  to  convert  the  x-ray  quanta  into  visible 
light  photons  that  are  subsequently  collected  by  an  optical  detector.  Characteristic  x  rays 
can  be  produced  within  the  phosphor.  Some  of  these  characteristic  x-ray  quanta  are 
reabsorbed  in  the  phosphor.  This  characteristic  x-ray  re-absorption  affects  the  quality  of 
images.  In  this  paper,  we  developed  a  3D  convolution  model  to  numerically  simulate  the 
characteristic  x-ray  re-absorption  inside  an  image  detector  or  intensifying  screen  with 
Monte  Carlo  technology.  The  3D  model  is  built  through  analyzing  the  physic  process  of 
creation,  traveling,  and  absorption  of  characteristic  x-ray  quanta  so  that  it  can  simulate 
the  noise  properties  of  characteristic  x-ray  re-absorption  accurately.  The  phosphor  layer 
of  an  image  detector  is  conceptually  divided  into  multiple  sub-layers  and  every  sub-layer 
is  conceptually  divided  into  multiple  voxels  based  on  the  pixel  size  of  a  digital  image. 

The  3D  model  computes  the  number  distribution  of  reabsorbed  characteristic  x-ray 
quanta  in  every  sub-layer  using  a  pre-calculated  convolution  kernel.  The  convolution 
kernel  gives  the  probability  distribution  that  a  characteristic  x-ray  quantum  can  be 
reabsorbed  by  a  voxel.  A  cut-off  window  of  convolution  is  calculated  to  reduce  the 
computation  time.  This  3D  convolution  model  shows  that  most  characteristic  x-ray  re¬ 
absorption  blurs  the  imaging.  However,  some  of  characteristic  x-ray  re-absorption 
enhances  the  primary  x-ray  input,  depending  on  the  relative  position  between  the  creation 
site  and  re-absorption  site  of  the  characteristic  x  rays.  The  simulation  result  is  used  to 
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quantitatively  analyze  the  effect  of  the  characteristic  x-ray  re-absorption  on  the  imaging 
performance  of  the  image  detector  or  intensifying  screen.  The  calculation  to  a  Cesium 
Iodide  layer  with  pixel  size  17x17  fim 2  shows  that  the  characteristic  re-absorption  may 
reduce  the  modulation  transfer  function  of  the  image  detector  as  much  as  10%  at  some 
spatial  frequency,  which  is  close  to  the  measurement.  The  model  of  the  characteristic  x- 
ray  re-absorption  in  this  paper  can  be  combined  to  the  model  that  simulates  the  whole 
imaging  process  of  an  image  detector  or  intensifying  screen  to  generate  computer 
simulated  x-ray  images. 

Keywords:  image  detector,  intensifying  screen,  characteristic  x  rays,  re-absorption, 
Monte  Carlo  simulation 

1.  Introduction 

If  the  energy  of  incident  x  rays  is  above  the  K-edge  of  the  phosphor  layer  inside  an 
image  detector  or  intensifying  screen,  characteristic  x-ray  radiation  can  be  produced. 
Generally,  these  characteristic  x  rays  are  emitted  isotropically  and  most  likely  reabsorbed 
at  the  projected  position  near  the  site  where  the  incident  x  rays  interact  with  the  phosphor 
layer.  However,  some  proportion  of  the  characteristic  radiation  is  reabsorbed  at 
substantial  distances  from  the  initial  interaction  sites  [1]  and  the  re-absorption  affects  the 
quality  of  the  images  [2-4],  Figure  1  schematically  shows  the  re-absorption  of 
characteristic  x  rays  inside  the  phosphor  layer  of  an  image  detector  or  intensifying  screen, 
which  blurs  images.  In  Figure  1,  characteristic  x  rays  are  generated  from  the  interaction 
of  incident  x  rays  with  the  phosphor  layer  in  the  voxel  centered  at  point  A.  These 
characteristic  x  rays  isotropically  emit  inside  the  phosphor  material.  Some  of  these 


characteristic  x  rays  are  reabsorbed  in  the  phosphor  voxel  centered  at  point  B,  which 
results  in  visible  light  photons  there.  This  means  that  Pixel  b  gets  some  visible  light 
output  because  of  the  incident  x  rays  at  pixel  a,  which  blurs  the  image. 

The  effect  of  characteristic  fluorescence  inside  the  phosphor  layer  of  an  image 
detector  or  intensifying  screen  on  images  has  been  researched  25  years  ago.  Vybomy  et 
al  emphasized  that  the  escape  of  characteristic  x  rays  from  the  intensifying  screens  should 
be  taken  into  account  when  the  efficiencies  of  screen-film  systems  were  determined  [5]. 
Arnold  and  Bjarngard  experimentally  observed  the  spread  of  line  spread  function  (LSF) 
caused  by  K-characteristic  x  rays  in  rare  earth  phosphor  intensifying  screens  [6],  Zhao  et 
al  found  that  the  effect  of  K-fluorescent  re-absorption  on  the  imaging  performance  of  a- 
Se  digital  mammography  detectors  can  be  significant  [7].  Chan  and  Doi  studied  the  x-ray 
energy  absorption  by  Monte  Carlo  simulation  [8].  There  are  also  some  investigators  who 
used  analytical  methods  to  qualitatively  calculate  the  energy  deposited  by  the  primary 
beam  and  by  the  reabsorbed  K-characteristic  x  rays  as  an  estimate  of  the  total  energy 
absorbed  in  the  fluorescent  intensifying  screens  [9-15].  All  these  investigations  took  the 
phosphor  layer  as  a  single  thin  layer  and  thus  assumed  the  characteristic  x  rays  transfer  in 
a  2  dimensional  plane.  In  this  paper,  we  propose  a  3D  convolution  model  to  calculate  the 
number  distribution  of  the  reabsorbed  characteristic  x  rays  and  a  method  to  simulate  the 
physical  process  of  characteristic  x-ray  re-absorption  with  Monte  Carlo  technology  in  3 
dimensions.  With  this  3D  convolution  model  and  simulation  method,  we  can  accurately 
simulate  the  noise  property  of  the  characteristic  re-absorption  and  quantitatively  calculate 
the  variation  in  modulation  transfer  function  of  an  image  detector  or  intensifying  screen 
caused  by  the  characteristic  x-ray  re-absorption.  This  3D  characteristic  x-ray  re- 


4 


absorption  model  can  be  combined  to  the  model  that  simulates  the  whole  imaging  process 
of  an  image  detector  or  intensifying  screen  to  generate  computer  simulated  x-ray  images. 

2.  The  model  to  calculate  the  number  of  the  reabsorbed  characteristic  x  rays 

There  are  3  assumptions  in  our  model  to  calculate  the  number  of  the  reabsorbed 
characteristic  x  rays  inside  the  phosphor  layer  of  an  image  detector  or  intensifying  screen. 

1)  The  phosphor  layer  of  an  image  detector  or  intensifying  screen  is  homogeneous. 

2)  The  emission  of  the  characteristic  x  rays  created  inside  the  phosphor  layer  is 
isotropic. 

3)  The  predominant  characteristic  x  rays  are  K-characteristic  x  rays.  The  characteristic 
x  rays  from  other  shells  such  as  L  shell  are  negligible  comparing  with  the  K- 
characteristic  x  rays. 

As  shown  by  Figure  1,  the  phosphor  layer  of  an  image  detector  or  intensifying  screen 
is  conceptually  divided  into  multiple  thin  sub-layers.  All  these  sub-layers  have  the  same 
thickness  A  t .  If  the  input  x-ray  fluence  to  the  input  interface  of  the  phosphor  layer 
is  S(x,  y,0) ,  the  x-ray  fluence  to  a  sub-layer  at  depth  t  will  be  5(x,  y,  t )  and 

S(x,y,t)=  S(x,y,0)e^'  (1) 

where  //  is  the  linear  attenuation  coefficient  of  the  incident  x  rays  inside  the  phosphor 
layer. 

If  the  pixel  area  of  an  image  is  dxdy  ,  the  number  of  x-ray  quanta  removed  from  the 
incident  beam  at  a  voxel  dxdydt  centered  at  (x,y,t)  is 


dN(x,y,t)  =  juS(x,y,t) dxdydt  =  jUS(x,y, 0)e  f" dxdydt 


(2) 


With  Equation  (2),  the  number  of  K-characteristic  x-ray  quanta  created  in  the  voxel 
centered  at  ( x ,  y,t )  can  be  expressed  as 


dnk  ( x ,  y,t )  =  0)dN(x ,  y,t )  =  %(dS(x,  y)pe~fU  dxdydt  (3) 

where  dnk  (x,  y,t)  is  the  number  of  K-characteristic  x-ray  quanta  created  at  phosphor 
voxel  centered  at  (x,  y,t)',£  is  the  probability  that  an  interacting  x  ray  undergoes  a  K- 
shell  interaction;  (0  is  the  fluorescent  yield  of  K-shell  photoelectric  interactions  and 
S(x,  y)  =  S(x,  y,0) . 

Since  the  K-characteristic  x  rays  created  inside  the  phosphor  layer  emit  isotropically 
from  the  sites  where  they  are  created,  according  to  Equation  (3),  the  fluence  of  the 
characteristic  x-ray  quanta  created  in  the  voxel  centered  at  (x,  y,  t)  is,  as  shown  by 

Figure  2, 

dnk(*,y.t)e-’”  (4) 

4  np 

The  number  of  the  characteristic  x-ray  quanta  created  in  the  voxel  centered  at  (x,  y,  t)  that 
pass  through  spherical  shell  1  in  the  phosphor  layer,  which  is  centered  at(x,y,t)and  with 
radius  p ,  will  be 

dn^ ,  .  d  ^  yJ)e-„f  (5) 

In  the  same  way,  we  can  get  the  number  of  the  characteristic  x-ray  quanta  created  at  (x,  y, 
t)  that  pass  through  spherical  shell  2,  which  is  with  radius  p  +  dp . 

nk  \p+dp  =  dnk  (x,  y,  (p+dp)  (6) 
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From  Equation  (5)  and  (6),  we  can  get  the  number  of  the  characteristic  x-ray  quanta 
created  in  the  voxel  centered  at  (x,  y,  t)  that  are  removed  from  the  fluence  in  the  volume 
between  spherical  shell  1  and  spherical  shell  2.  This  number  is 

Ank  =  nk  -nk\p+dp=  dnk (x, y,t){e~flkP  —  e~f‘i(p+dp) }  (7) 

In  Equation  (7),  if  dp  — >  0 ,  the  density  of  the  characteristic  x-ray  quanta  created  in  the 

voxel  centered  at  (x,  y,  t)  in  the  volume  between  spherical  shelll  and  spherical  shell  2  is 

A nk  _  dnk(x,y,t)(e'/JkP  - e~f‘tip+dp) )  _  dnk(x,y,t)  d  ^ 

AV  Anp2dp  Anp 2  dp 


According  to  Equation  (8),  we  can  get  the  probability  that  the  voxel  centered  at  (x,  y,t) 
can  reabsorb  a  characteristic  x-ray  quantum  created  in  the  voxel  centered  at  (x,  y,  t),  as 
shown  by  Figure  3, 


p(x',y',t',x,  y,t )  = 


1 


dnk(x,y,t) 


dnk(x,y,t )  4  np 


1 


—  (e'MtP)dv'^— —pke-^pdv'  (9) 

dp  Anp1  k 


where  rfv'is  the  volume  of  the  voxel  centered  at(x',  y',t') ;  p  is  the  distance  between  the 
creation  phosphor  voxel  of  the  characteristic  x-ray  quanta  and  the  observed  phosphor 
voxel  that  reabsorbs  some  of  these  characteristic  x-ray  quanta  and 

p=J(x'-x)2  +  (y'-y?  +  (t'-t)2  (10) 


If  /?  — >  0,  it  means  that  the  characteristic  rays  travel  an  infinitesimal  distance  and 
attenuation  is  zero,  thus  p{x',  y',t',x,  y,t)  — >  0 because  there  is  no  characteristic  quantum 
removed  from  the  characteristic  fluence  in  zero  distance.  If  we  use  dmk  (x\  y',t',x,  y,t )  to 
denote  the  number  of  the  characteristic  x-ray  quanta  created  in  the  voxel  centered  at  (x,  y, 
t)  and  reabsorbed  by  the  voxel  centered  at(x',  y',t'),  using  Equation  (9),  this  number  can 


be  expressed  as: 
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dmk(x',y',t',x,y,t)  =  p(x',y',t',x,y,t)dnk(x,y,t)  (11) 

In  fact,  every  voxel  in  the  phosphor  layer  of  an  image  detector  or  intensifying  screen  may 
produce  some  characteristic  x-ray  quanta.  Using  Equation  (3)  and  Equation  (11),  the  total 
number  of  reabsorbed  characteristic  x-rays  by  the  voxel  centered  at  (x\  y',t')  can  be 
expressed  as 

mk  (x',y',t')  =  ^p(x\y\t\x,y,tWx,yfi)c*dxdydt  (12) 

wholescreen 

where  mk  (x\  y',t')  is  the  total  number  of  reabsorbed  characteristic  x-rays  by  the  voxel 
centered  at  (x',  y',t') .  The  integration  in  Equation  (12)  is  over  the  whole  phosphor  layer. 

To  calculate  the  integration  of  Equation  (12)  numerically,  we  can  conceptually 
divide  the  phosphor  layer  into  multiple  sub-layers  that  have  the  same  thickness  At .  This 
way,  the  whole  phosphor  layer  becomes  conceptually  a  set  of  voxels  that  have  equal 
volume  AxAyAt ,  where  AxAy  is  the  area  of  an  image  pixel.  Then,  the  integration  of 
Equation  (12)  can  be  written  as 

mk(i',j',k')  =  j)p(i',  j',k',i,  AxAyAt  (13) 

i  j  * 

where  ( i ,  j,k)  is  the  index  of  a  source  vexel  where  the  characteristic  x-ray  quanta  are 
created;  (/',  j\k')  is  the  index  of  the  observed  voxel  that  reabsorbs  the  characteristic  x-ray 
quanta  from  all  other  voxels;  S(i,  j) is  the  input  x-ray  quantum  distribution  to  the  input 
interface  of  the  phosphor  layer  of  an  image  detector  or  intensifying  screens  and 
S(i,  j)  =  S(x,y)S(iAx,  jAy) ,  where  S(iAx,jAy) is  a  2D  delta  function  and  S(x,y)  is  the 
primary  x-ray  input  to  image  pixel  (i,  j ) ;  p{i',  j',k',i ,  j,k)  is  the  probability  that  a 
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characteristic  x-ray  quanta  created  in  voxel  O',  j, k)  is  reabsorbed  by  voxel  ( i',j',k ');  Nx, 
N  ,and  Nt  are  total  voxel  numbers  in  x,  y,  and  t  direction. 

3.  The  reabsorbed  characteristic  x-rays  that  blur  the  image  and  enhance  the 
primary  x-ray  input 

In  order  to  analyze  the  effect  of  the  reabsorbed  characteristic  x-ray  quanta  on  the 
imaging  system,  we  can  separate  the  integration  domain  of  equation  (13)  into  2  parts,  as 
shown  by  Figure  4.  To  any  observed  voxel  (i\j\  k') ,  domain  A  includes  the  voxels  with 
i  ^  i'  or  j  ^  j'  and  domain  B  includes  the  voxels  with  i  =  i'  and  j  =  / .  In  fact,  domain  B 
is  a  phosphor  bar  from  pixel  (/',  j', 0)  on  the  input  interface  to  (/',  j',T)  on  the  bottom  of 
the  phosphor  layer,  and  domain  A  includes  all  other  voxels  of  the  phosphor  layer  except 
the  voxels  of  that  phosphor  bar. 

The  characteristic  x-ray  quanta  created  in  the  voxels  of  domain  A  are  not  related  to 
the  primary  x-ray  input  to  pixel  (T,  j', 0)  and  do  not  contain  image  information  to  observed 
voxel  (/',  /,  k’)  so  that  these  characteristic  x-ray  quanta  blur  pixel  (i\  j',0 ) .  Then,  we  can 
call  the  re-absorption  of  these  characteristic  x-ray  quanta  as  blurring  characteristic  x-ray 
re-absorption  or  crossing  characteristic  x-ray  re-absorption  because  these  characteristic  x- 
ray  quanta  must  have  crossed  at  least  one  pixel  boundary  before  arriving  at 
voxel  (/',  j',k') .  Re-absorption  of  characteristic  x-ray  quanta  created  in  domain  A  is 
equivalent  to  voxel  (/',  j',k')  s’  obtaining  some  x-ray  quanta  from  many  blurring 


background  x-ray  sources. 


9 


The  characteristic  x-ray  quanta  created  in  the  voxels  of  domain  B  are  due  to  the 
primary  x-ray  input  to  pixel  (/',  j',0 ) .  These  characteristic  x-ray  quanta  do  not  cross  any 
pixel  boundary  when  arriving  at  the  observed  voxel  (T,  j\  k')  ,  as  shown  by  Figure  5.  The 
primary  x-ray  input  at  pixel  (i',  j',0)  contains  both  image  information  and  noise.  But,  the 
x-ray  quanta  containing  image  information  and  the  x-ray  quanta  containing  noise  have 
the  same  energy  distribution  and  have  the  same  probability  to  create  characteristic  x  rays. 
If  there  are  more  primary  input  x-ray  quanta  containing  the  image  information,  the  re¬ 
absorption  of  the  characteristic  x-ray  quanta  created  in  domain  B  enhances  the  image 
information.  If  there  are  more  noise  quanta  in  the  primary  x-ray  input  to  pixel  (/',  j',0) , 
the  re-absorption  of  the  characteristic  x-ray  quanta  created  in  domain  B  enhances  noise. 
However,  no  matter  if  the  characteristic  x-ray  quanta  are  created  from  noise  x-ray  input, 
re-absorption  of  characteristic  x-ray  quanta  created  in  domain  B  is  equivalent  to  voxel 
( i ',  j',k ')  ’s  getting  more  x-ray  quanta  from  primary  x-ray  input  at  pixel  ii',  j',0)  so  that 
we  can  call  this  characteristic  x-ray  re-absorption  as  enhancing  characteristic  x-ray  re¬ 
absorption. 

Based  on  the  analysis  above,  we  can  see  that  there  are  three  types  of  x-ray  input 
quanta  to  every  voxel  of  phosphor  layer  when  including  characteristic  x-ray  re-absorption: 
the  absorbed  primary  input  x-ray  quanta,  the  reabsorbed  characteristic  x-ray  quanta 
enhancing  the  primary  x-ray  input,  and  the  reabsorbed  characteristic  x-ray  quanta  adding 
the  blurring  background. 

4.  Convolution  model  of  characteristic  x-ray  re-absorption 


From  Equation  (9),  we  can  have 
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p(i',  j',k\ i,  j,  k)  =  - — - - 

4 


p(i  ,j,k  ,i,j,k)  =  0  when  pr,yx,u,k  =  0 
Ax'  =  Ax ,  Ay'  =  Ay  ,  At'  =  At 


-MkPr.fyj.jA  Av'Av'A/ 


Ajf  Ay  At  when  Av,*',/.;,*  *  0 


(15) 


(14) 


Pi\j\k\i.j.k  =  V(A-02  {i'-if  +  (Ay)2  (/  -  jf  +  (At)2  {k'-kf  (16) 


Using  Equation  (14)  to  (16),  Equation  (13)  can  be  written  into 


mk(i',j',k')  = 
^0)(AxAyAtf  ppk  ^ 


4  n 


lllsa.j)1 


<•  j  A 


-(tikAi+jJtPr.i-x.i.j.k ) 


Pi'.j’.k'.i.j.k 


Pi\rx.ij.k  *  0 


(17) 


mk (i ,  j  ,k  )  =  0 


/i',j’,k',i,j,k 


=  0  (18) 


Since  S(i,  j)  is  not  a  function  of  index  k,  we  can  rewrite  Equation  (17)  into 


. .  _  fdAwufm  *•&...  ^ 


mk (i ,  j  ,k)  = 


4  n 


■11  S(i.j)l 


k  Pi',j',k',i,j,k  J 


(19) 


From  Equation  (16),  we  know  that  Pyyx,ij,k  Is  or|ly  a  function  of  (/'  -  i,  /  -  &) . 

Then,  the  number  distribution  of  reabsorbed  characteristic  x-ray  quanta  by  a  sub-layer  of 
phosphor  with  index  k'can  be  expressed  as  a  convolution  as  shown  by  Equation  (20). 

Nx  Ny 

mk  (k')  =  X  2 Sd,  j)WbUir  a'  -  i,  j'  -j,k')  =  s*  W(k')  (20) 

i  j 


where  mk  ( k ')  is  the  number  distribution  of  reabsorbed  K-characteristic  x-ray  quanta  in 
sub-layer  k' ,  W(i'-i,  j' -  j,k') is  the  convolution  kernel  and  it  is  the  characteristic  x-ray 
re-absorption  probability  function  and  defined  as 
Za){AxAyAt)2 ppk  A 

W(i  -i,j -j,k)  =  *  '  ~  Z - 5 -  when  Pwjj*  *° 

^  A  Pi'.y.k’.i.j.k 
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w  =  0  when  pftfXMJc  =  0 


or 


........  <^O^AxAyAt)2  HMk  ^  e~fAA‘~Mii  .... 

W(i,j,k')  =  — - —  £ - 3 -  when  Pww*  0  (21) 


An 


T  Ptj,k 


W(  0, 0,0)  =  0 

From  Section  3,  we  know  that  the  reabsorbed  characteristic  x-ray  quanta  have 
different  imaging  performance  depending  on  the  relative  position  of  the  re-absorption  site 
to  the  creation  site  of  the  reabsorbed  characteristic  x-ray  quanta.  The  elements  of 
convolution  kernel  W(i,  j,k ')  with  i  *  0  or  j  *  0  are  related  to  the  characteristic  x-ray 
quanta  created  in  the  domain  A  of  Figure  4.  We  can  define  these  elements  as  a 
characteristic  blurring  probability  function. 

Wb,lirri„g(k')  =  W(i,j,k')  when  i*0orj*0  (22) 

The  elements  of  convolution  kernel  W (0,0,  k')  are  related  to  the  characteristic  x-ray 
quanta  created  in  domain  B  of  Figure  4;  we  can  define  these  elements  as  characteristic 
enhancing  probability  factors. 

Wenhancing(k')  =  W(0,0,k')  (23) 


Substitute  i  with  0  and  j  with  0  in  Equation  (16)  and  Equation  (21),  we  can  get 


till  N‘  p-VkM  -vk\k'-k\M 


(24) 


With  Equation  (22)  and  Equation  (23),  we  have 

^(^')=ww„m,s(r)+w£nW„s(fc') 


(25) 


According  to  Equation  (21),  the  characteristic  x-ray  re-absorption  probability  function 
is  decided  by  the  linear  attenuation  coefficient  of  the  primary  input  x  rays,  the  linear 
attenuation  coefficient  of  the  characteristic  x  rays,  the  spatial  resolution  of  digital  images, 
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and  the  number  of  sub-layers  that  the  whole  phosphor  layer  is  conceptually  divided.  To  a 
given  image  detector  or  intensifying  screen,  the  characteristic  x-ray  re-absorption 
probability  function  of  its  phosphor  layer  can  be  calculated  using  Equation  (25)  for  every 
sub-layer  whenever  the  pixel  size  and  sub-layer  thickness  are  confirmed. 

Theoretically,  the  convolution  of  Equation  (20)  should  overlap  the  whole  phosphor 
voxel  set  and  the  element  number  of  characteristic  x-ray  re-absorption  probability 
function  of  every  sub-layer  should  be  as  same  as  the  pixel  number  of  the  image.  This  will 
require  not  only  a  lot  of  time  to  compute  the  characteristic  x-ray  re-absorption  probability 
function,  but  also  much  more  time  to  do  the  convolution  for  every  sub-layer  when 
calculating  the  number  distribution  of  characteristic  x-ray  quanta  reabsorbed,  especially  if 
the  image  has  many  pixels  and  the  phosphor  layer  of  the  image  detectors  or  intensifying 
screen  is  conceptually  divided  into  many  very  thin  sub-layers.  However,  as  characteristic 
x-ray  quanta  are  continuously  removed  from  the  beams  when  traveling  inside  the 
phosphor  layer,  the  fluence  intensity  of  the  characteristic  x  rays  created  inside  the 
phosphor  layer  decays  exponentially.  Then,  the  probability  that  the  characteristic  x-ray 
quanta  are  reabsorbed  by  the  voxels  relatively  far  from  the  creating  sites  drops 
exponentially.  Table  1  lists  the  distances  that  the  K-characteristic  x  rays  can  travel  inside 
the  phosphor  layer  before  the  intensity  of  fluence  is  attenuated  to  1%,  5%,  and  10%  of  the 
original  intensity  inside  Gadolinium  Oxysulfide  and  Cesium  Iodide.  The  distance  values 
in  Table  1  are  calculated  using  the  attenuation  coefficient  of  the  x  rays  at  K-edge 
energies  of  the  phosphor  materials.  If  we  ignore  the  minor  contribution  of  the 
characteristic  x  rays  that  are  created  relatively  far  from  the  observed  voxel,  the 
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convolution  of  Equation  (20)  only  needs  cover  a  limited  number  of  voxels  around  the 
observed  voxel  as  shown  by  Figure  6. 

To  any  voxel  centered  at  the  bar  from  ( x ',  y',0)  to  (x\  y',  T)  in  Figure  6,  the 
computation  to  calculate  the  number  of  characteristic  x  ray  quanta  needs  only  to  cover  a 
cubic  2n  pixels  long,  2n  pixels  wide,  and  as  thick  as  the  phosphor  layer  shown  by  the 
window  in  Figure  6.  We  call  this  window  as  a  cut-off  window.  The  n  of  a  cut-off 
window  can  be  decided  according  to  the  phosphor  materials  used  in  the  image  detector  or 
intensifying  screen,  the  size  of  the  image  pixels,  and  the  expected  accuracy  of  the 
calculation.  Table  2  lists  the  pixel  numbers  that  K-characteristic  x  rays  can  travel  inside 
the  phosphor  materials  before  99%,  95%,  or  90%  quanta  are  removed  from  the  beams  at 
pixel  size  50x50//m2and  17 xlljUm2 for  Gadolinium  Oxysulfide  and  Cesium  Iodide. 
With  this  table,  we  can  choose  the  size  of  a  cut-off  window.  For  example,  if  the 


Table  1  Physical  characteristic  of  two  phosphor  materials 


Phosphor 

Materials 

K-edge 

(KeV) 

Density 
(g  /cm3) 

Attenuation 
(cm2/ g)  of  K 
characteristic 
x  rays 

Distance 
to  1%  of 
the 

original 

(cm) 

Distance 
to  5%  of 
the 

original 

(cm) 

Distance 
to  10% 
of  the 
original 
(cm) 

Gd202S 

50.2 

7.34 

3.234 

0.128 

0.098 

Csl 

36.0/33.2 

4.51 

6.923 

0.147 

0.096 

0.074 

calculation  wants  to  include  95%  characteristic  x-ray  quanta  created  inside  the  phosphor 
layer  of  an  image  detector  or  intensifying  screen,  the  cut-off  window  needs  to  cover  55 
pixels  in  one  of  x  or  y  directions  when  using  Cesium  Iodide.  Then,  the  n  in  the  window  of 
Figure  5  can  be  set  as  55.  The  kernel  of  the  convolution  in  Equation  (20),  characteristic 
x-ray  re-absorption  probability  function,  needs  to  have  1 1 1  elements  in  x  dimension  and 
y  dimension,  with  1 10  elements  of  Wblurring  and  1  element  of  Wenhancing . 
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Table  2  Cut-window  size  to  calculate  the  convolution  of  Equation  (20) 


Phosphor 

Materials 

Pixel  size  50x50 /Jm2 

Pixel  size  17xl7//m2 

Including 
99%  of 
Characte¬ 
ristic  X 
quanta 

Including 
95%  of 
Characte¬ 
ristic  X 
quanta 

Including 

90%  of 
Characte¬ 
ristic  X 
quanta 

Including 

99%  of 
Characte¬ 
ristic  X 
quanta 

Including 

95%  of 
Characte¬ 
ristic  X 
quanta 

Including 
90%  of 
Characte¬ 
ristic  X 
quanta 

Gd202S 

40 

26 

20 

113 

74 

56 

Csl 

30 

20 

15 

84 

55 

43 

5  Quantitatively  evaluate  the  effect  of  the  reabsorbed  characteristic  x  rays  on 
imaging  performance 

Using  Equation  (20)  and  Equation  (21),  we  can  simulate  the  characteristic  x-ray  re¬ 
absorption  process  inside  a  phosphor  layer  and  quantitatively  evaluate  the  affect  of 
characteristic  x-ray  re-absorption  to  the  image  performance  of  an  image  detector. 

a)  Calculate  the  number  of  reabsorbed  characteristic  x-ray  quanta  blurring  the 
image. 

We  know  that  some  reabsorbed  characteristic  x-ray  quanta  do  not  contain  any  image 
information  and  only  blur  the  image  and  some  reabsorbed  characteristic  x-ray  quanta 
carry  image  information  and  enhance  the  primary  input.  A  necessary  step  to  calculate  the 
number  of  reabsorbed  characteristic  x-ray  quanta  blurring  the  image  by  every  phosphor 
voxel  is  to  calculate  the  convolution  kernel  of  Equation  (22),  the  characteristic  x-ray 
blurring  probability  function  Wblurring  (k')  for  every  sub-layer  of  the  phosphor  material 
according  to  the  pixel  size  and  the  selected  accuracy.  Using  Equation  (21)  and  Equation 
(22),  Wblurring  (/,  j,k')  can  be  calculated  as  a  three  dimensional  table.  Figure  7  and  Figure 

8  are  two  plots  of  the  characteristic  x-ray  blurring  probability  function  Wblur  (/,  j,k')  for 
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the  first  sub-layer  of  phosphor  (the  sub-layer  is  mostly  close  to  the  incident  interface  of 
the  phosphor  layer.)  These  two  plots  are  the  calculation  results  for  a  Cesium  Iodide 
phosphor  layer  [17].  The  total  thickness  of  the  Cesium  Iodide  layer  is  assumed  as  85  /urn. 
When  doing  the  calculation,  the  whole  Cesium  Iodine  layer  is  conceptually  divided  into 
40  sub-layers  and  all  sub-layers  have  the  same  thickness,  2.125  Jilin .  The  image  pixel  size 

is  assumed  asl7xl7//m2so  that  the  size  of  every  Cesium  Iodide  voxel  isdx  =  \l  Jim, 
dy  =  17  jLUn ,  and<*  =  2.125  jJm.  The  accuracy  of  calculation  is  set  95%,  which  includes 
95%  characteristic  x-ray  quanta  created  in  calculation.  Then,  according  to  Table  2,  the 
size  of  the  cut-off  window  in  Figure  6  is  55  and  Whlur  (/,  j,  k)  is  a  40x1 10x1 10  table.  The 
energy  of  incident  x  rays  is  assumed  as  53  keV.  At  this  energy,  the  linear  attenuation 
coefficient  of  Cesium  Iodide  is  fJ.  =  51.34  cm-1 .  The  K-edge  energy  of  Cesium  Iodide  is 
36.0/33.2  keV  [1].  The  attenuation  coefficient  at  this  energy  is  juk  =  32.2  cm' .  The 
probability  that  an  interacting  x  ray  undergoes  a  K-shell  interaction,  is  set  as  0.85  and 
the  fluorescent  yield  of  K-shell  photoelectric  interactions,  Q) ,  is  set  as  0.93  according  to 
[13].  Figure  7  is  the  contour  of  the  characteristic  blurring  probability  function  of  the  first 
sub-layer.  Figure  8  is  the  relative  distribution  of  the  characteristic  blurring  probability 
function  of  the  first  sub-layer.  The  two  plots  show  that  most  of  the  reabsorbed 
characteristic  x-ray  quanta  blurring  the  image  are  from  the  voxels  that  are  close  to  the 
observed  voxel.  This  is  because  the  inverse  square  law  shown  by  Equation  (21).  The 
intensity  of  the  characteristic  x  rays  becomes  smaller  when  the  observed  voxel  is  farther 
from  the  creation  site  of  the  characteristic  quanta.  Equation  (21)  also  shows  the 
probability  of  re-absorption  is  proportional  to  the  square  of  the  voxel  volume.  Since  the 
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voxel  volume  in  our  calculation  is  only  614.125  /zm3 ,  the  absolute  values  of  the 

characteristic  x-ray  blurring  probability  function,  which  have  unit  mm 2 ,  are  relatively 
small  as  shown  by  Figure  7. 

b)  Calculate  the  number  of  reabsorbed  characteristic  x-ray  quanta  enhancing  the 
primary  x-ray  input. 

Figure  9  is  the  plot  of  the  characteristic  x-ray  enhancing  probability  factor  vs.  the 
sub-layer  index.  These  characteristic  x-ray  enhancing  probability  factors  are  calculated 
with  Equation  (24)  using  the  same  parameters  as  used  when  calculating  the  characteristic 
x-ray  blurring  probability  functions.  In  Figure  9,  the  deeper  the  sub-layer  is,  the  bigger 
the  index  is.  The  sub-layer  with  the  input  interface  of  the  phosphor  layer  has  an  index  0 
and  the  bottom  sub-layer  has  an  index  39.  The  value  of  characteristic  x-ray  enhancing 
probability  factor  is  not  very  big  because  this  value  is  proportional  to  the  square  of  pixel 
area. 

From  Figure  9,  we  can  see  that  there  is  a  build-up  stage  in  characteristic  x-ray 
enhancing  probability  factor.  The  first  sub-layer  has  smaller  characteristic  x-ray 
enhancing  probability  factor.  It  increases  with  the  sub-layer  index  (deeper  in  the 
phosphor  layer).  After  arriving  at  a  peak  value,  the  characteristic  x-ray  enhancing 
probability  factor  decreases  with  sub-layer  index.  This  is  because  the  first  sub-layer  can 
only  reabsorb  characteristic  x-ray  quanta  created  beneath  it.  However,  the  sub-layer 
deeper  in  the  phosphor  layer  can  reabsorb  the  characteristic  x-ray  quanta  created  both 
above  it  and  below  it,  which  makes  the  characteristic  x-ray  enhancing  probability  factor 
increases.  Because  of  the  exponential  factor  in  Equation  (24),  the  characteristic  x-ray 
factor  decreases  as  the  sub-layer  is  far  from  the  input  interface  of  the  phosphor  layer.  So, 
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after  the  peak  point,  as  shown  by  Figure  9,  the  characteristic  x-ray  enhancing  probability 
factor  decreases  with  depth.  The  characteristic  x-ray  factor  drops  quickly  when  the  sub¬ 
layers  are  close  to  the  bottom  of  the  phosphor  layer  since  the  voxels  on  these  sub-layers 
can  only  reabsorb  the  characteristic  x-ray  quanta  created  above  them, 
c)  Calculate  the  variation  in  modulation  transfer  function  caused  by  characteristic 
x-ray  re-absorption  by  numerical  simulation 

One  way  to  research  the  effect  of  the  characteristic  x-ray  re-absorption  on  the  imaging 
performance  is  to  analyze  the  influence  of  the  characteristic  x-ray  re-absorption  to  the 
modulation  transfer  function  of  an  image  detectors  or  intensifying  screen  [18-22],  Using 
the  calculated  characteristic  x-ray  blurring  probability  function  and  enhancing  probability 
factor,  we  can  simulate  the  characteristic  x-ray  re-absorption  process  with  Monte  Carlo 
technique.  With  the  simulation  output,  we  can  calculate  the  modulation  transfer  function 
of  the  image  detector  or  intensifying  screen.  We  did  this  simulation  in  following  steps 
referring  to  Cunningham’s  cascaded  linear  system  analysis  approach  [23-25]. 

Step  1  Building  an  x-ray  image  input. 

As  we  know,  to  obtain  the  modulation  transfer  function  of  an  image  detector,  we 
need  a  point  x-ray  image  input  to  the  image  detecting  system.  We  used  Equation  (26)  as 
our  image  input.  This  is  an  ideal  point  input  because  there  is  no  x-ray  quantum  to  all 
other  pixels  except  the  central  pixel.  The  x-ray  input  is  assumed  monochromic  and  at 
energy  53  keV. 


1525  x  -  ray  quanta  when  i  =  127,  and  j  =  127 
S(i';',  =  l  0  else  (26) 


where  i  and  j  are  pixel  indices  of  the  image.  We  simulated  an  image  detector  that  has 


256  pixels  in  width  and  256  pixels  in  height.  The  pixel  size  of  the  image  is 
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17x17  // m2  so  that  the  input  fluence  at  the  center  of  the  image  has  a  pulse  of5.3xl06 
quanta!  mm2 

Step  2  Calculating  the  numbers  of  reabsorbed  characteristic  x-ray  quanta  blurring  the 
image  and  enhancing  the  primary  input. 

This  calculation  can  be  done  for  every  pixel  using  the  convolution  model,  Equation 
(20)  to  (24),  with  the  pre-calculated  character  x-ray  blurring  probability  function  and 
enhancing  probability  factor.  Through  this  calculation,  we  obtained  the  number 
distributions  of  reabsorbed  characteristic  x-ray  quanta  blurring  the  image  and  enhancing 
the  primary  x-ray  input  for  every  sub-layer.  Figure  10  plots  the  accumulated  number 
distribution  of  reabsorbed  characteristic  x-ray  quanta  blurring  the  image  under  the  ideal 
point  input  in  order  to  show  that  the  pixels  close  to  the  central  pixel  are  blurred  by 
reabsorbing  the  characteristic  x-ray  quanta  created  in  the  voxels  bellow  the  central  pixel. 
In  fact,  the  re-absorption  takes  place  in  every  sub-layer  so  that  the  re-absorbed 
characteristic  x-ray  quanta  have  different  effect  on  the  final  output  image  because  the 
visible  light  photons  produced  from  these  reabsorbed  characteristic  x-ray  quanta  have 
different  diffusion  and  absorption.  In  the  following  steps,  we  simulate  the  visible  light 
photon  production,  absorption,  and  diffusion  sub-layer  by  sub-layer. 

Table  3  lists  the  accumulated  number  of  the  reabsorbed  characteristic  x-ray  quanta 
blurring  the  image  for  the  image  pixels  close  to  the  image  center  where  the  primary  input 
x-ray  fluence  pulse  is.  One  problem  is  that  the  Monte  Carlo  simulation  outputs  are  real 
numbers  as  shown  in  Table  3.  Here,  the  decimal  part  can  be  explained  as  the  probability 
of  another  characteristic  x-ray  quantum  might  be  absorbed  by  the  observed  voxel.  For 
example,  the  pixel  near  the  center  pixel  has  a  number  3.668365.  This  can  be  understood 
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as  that  the  phosphor  voxels  bellow  this  pixel  reabsorbs  3  characteristic  x-ray  quanta. 
Besides,  there  might  be  the  fourth  characteristic  x-ray  quantum  reabsorbed  there.  The 
possibility  that  the  voxels  bellow  this  pixel  reabsorbs  the  fourth  characteristic  x-ray 
quantum  is  0.668365. 

From  Table  3,  we  can  see  that  the  pixels  close  to  central  pixel  will  have  some  image 
output  though  there  is  no  primary  x-ray  input  to  these  pixels.  The  characteristic  x-ray  re¬ 
absorption  is  equivalent  to  extra  x-ray  input  sources  to  the  phosphor  layer  of  an  image 
detector  or  intensifying  screen.  After  adding  reabsorbed  characteristic  x-ray  quanta  from 
the  Monte  Carlo  simulation,  every  sub-layer  of  the  phosphor  layer  gets  three  sources  of 
x-ray  inputs,  the  primary  x-ray  input,  the  blurring  x-ray  input  that  consists  of  reabsorbed 
characteristic  x-ray  quanta  blurring  the  image,  and  the  enhancing  x-ray  input  that  consists 
of  reabsorbed  characteristic  x-ray  quanta  enhancing  the  primary  x-ray  input  at  central 
pixel.  All  of  these  three  types  of  x-ray  inputs  may  generate  visible  light  and  contribute  to 
the  image  output. 

Step  3,  Simulate  the  conversion  process  from  x-ray  quanta  to  visible  light  photons 
inside  the  phosphor  layers. 


Table  3  The  number  distribution  of  reabsorbed  characteristic  x-ray  quanta  blurring  the 
image  under  the  ideal  point  x-ray  input. _ _ _ _ 


0.344182 

0.470862 

0.600216 

0.600216 

0.470862 

0.344182 

0.470862 

0.732893 

tmwM 

■irilUliM 

1.094941 

0.732893 

0.470862 

0.600216 

1.094941 

3.668365 

2.240517 

1.094941 

0.600216 

0.660068 

1.313550 

3.668365 

0.000000 

3.668365 

1.313550 

0.660068 

0.600216 

1.094941 

2.240517 

3.668365 

2.240517 

1.094941 

0.600216 

0.470862 

0.732893 

1.094941 

1.313550 

1.094941 

0.732893 

0.470862 

0.344182 

0.470862 

0.600216 

0.660068 

0.600216 

0.470862 

0.344182 

The  numbers  listed  are  accumulated  reabsorbed  characteristic  x-ray  quanta  from  all  sub-layers.  The  center 
of  the  table  is  the  center  of  the  image  where  the  input  x-ray  pulse  is.  The  phosphor  material  is  Cesium 
Iodide.  The  point  input  is  1525  x-ray  quanta. 
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We  use  0.15  as  the  conversion  efficiency  and  2.4  eV  as  the  energy  of  a  visible  light 
photon.  Then,  the  conversion  from  x-ray  quanta  to  visible  light  photons  can  be  simulated 


nB  ~  0.15 — —n 

p  2.4  * 


where  np  is  the  number  of  visible  light  photons  converted  from  x-ray  quanta;  nx  is  the 
number  of  x-ray  quanta;  Ex  is  the  energy  of  an  x-ray  quantum.  In  our  simulation,  Ex  is  53 
keV  for  primary  x-ray  input  and  34.5  keV  (average  of  36/33  keV)  for  characteristic 
blurring  input  or  enhancing  input. 

Step  4  Simulate  the  diffusion  and  absorption  of  the  visible  light  photons  inside  the 
phosphor  layer. 

The  diffusion  and  absorption  of  the  visible  light  photons  inside  the  phosphor  layer 
can  be  simulated  by  Boltzmann  equation  [26-30].  We  used  the  solution  of  Boltzman 
equation  in  spatial  frequency  space,  as  shown  by  Equation  (28). 


?Vv(M’v’0  = 


rp^a+rp^e0*  +  (0-t/7o)<Tct] 

{o)+tp0  fco  +  rp ,  y*  -  (a>-  xp0  ){a>- zpx  )e 


where  T  is  the  inverse  relaxation  length  that  describes  the  probability  of  optical 
absorption  occurring;  zw  is  related  to  spatial  frequency  and  (O2  =  <7 2  +  \n{u 2  +v2)  where 
cris  the  reciprocal  diffusion  length  and  characterizes  optical  scattering;  u  is  the  spatial 
frequency  in  horizontal  direction  and  v  is  the  spatial  frequency  in  vertical  direction;  t  is 
the  depth  of  the  sub-layer;  T  is  the  thickness  of  the  phosphor  layer  of  an  image  detector  or 


intensifying  screen;  /?,.= - where  r;  is  the  reflectivity  of  the  boundary  of  the  screen; 

1  +  r 


( i  =  0 )  represents  the  interface  closest  to  the  x-ray  source  of  the  screen  and  (  i  =  1) 
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represents  the  output  interface.  To  use  Equation  (28),  we  first  did  a  Fourier 
transformation  to  the  visible  light  photon  distribution  of  a  sub-layer,  obtained  in  the 
conversion  process  from  x-ray  quanta  to  visible  light  photons.  Then,  we  used  this  Fourier 
transformation  to  multiply  (pu  V(w,v,  t )  of  the  Equation  (28).  Finally,  we  did  an  inverse 

Fourier  transformation  to  get  the  visible  light  image  output. 

Our  simulation  flow  is  shown  by  Figure  11.  We  first  pickup  a  sub-layer  and 
calculate  the  number  distribution  of  absorbed  primary  input  x-ray  quanta.  Based  on  the 
pre-calculated  characteristic  x-ray  blurring  probability  and  enhancing  probability,  we 
calculate  the  number  distribution  of  the  reabsorbed  characteristic  x-ray  quanta.  After  this 
computation,  we  simulate  the  visible  light  photon  production,  absorption,  and  diffusion  of 
the  sub-layer  and  get  the  image  output  of  this  sub-layer.  Then,  we  pickup  another  sub¬ 
layer  and  do  the  same  calculation  and  simulation  again  until  we  get  the  image  output  of 
every  sub-layer.  Finally,  we  add  the  image  output  of  every  sub-layer  together  to  obtain 
the  image  output  of  the  whole  phosphor  layer  of  the  image  detector  or  intensifying  screen. 

By  the  simulation  processes  above,  we  obtained  a  simulated  image  output  of  the  ideal 
point  x-ray  input  image  from  a  Cesium  Iodide  layer.  According  to  the  definition  of 
modulation  transfer  function,  we  can  get  the  modulation  transfer  function  of  the  Cesium 
Iodide  layer  by  doing  Fourier  Transformation  to  the  simulated  image  output  and 
normalizing  the  Fourier  Transformation  output.  Figure  12  shows  the  modulation  transfer 
functions  calculated  for  the  Cesium  Iodide  layer  with  thickness  85  /tin  and  pixel 

sizel7xl7//m2 .  The  dashed  curve  in  Figure  12  is  the  modulation  transfer  function 
including  the  characteristic  x-ray  re-absorption.  The  solid  curve  is  the  modulation  transfer 
function  ignoring  the  characteristic  x-ray  re-absorption.  From  this  comparison,  we  can 
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see  that  the  spatial  resolution  of  the  image  detector  or  intensifying  screen  is  reduced  by 
the  characteristic  x-ray  re-absorption.  The  value  of  modulation  transfer  function  is 
reduced  from  spatial  frequency  0.1  mm~l  to  8  mm-1 .  The  biggest  reduction  is  about  10%, 
which  is  same  as  the  measured  variation  of  modulation  transfer  function  caused  by 
characteristic  x-ray  re-absorption  [6]. 

6.  Summary 

A  convolution  model  is  built  to  simulate  the  characteristic  x-ray  re-absorption  based 
on  analyzing  the  physics  process  of  characteristic  x-ray  production  and  transfer  inside  a 
phosphor  layer.  The  convolution  kernel  can  be  pre-calculated  based  on  the  physical 
properties  of  the  phosphor  material,  the  image  pixel  size,  and  the  selected  computation 
accuracy.  According  to  the  effect  on  the  imaging  performance,  the  reabsorbed 
characteristic  x-ray  quanta  inside  a  phosphor  layer  can  be  separated  into  two  types:  the 
reabsorbed  characteristic  x-ray  quanta  blurring  images  and  the  reabsorbed  characteristic 
x-ray  quanta  enhancing  primary  x-ray  input.  The  convolution  kernel  can  be  expressed  as 
a  summation  of  the  characteristic  x-ray  blurring  probability  function  and  the 
characteristic  x-ray  enhancing  probability  factor.  Using  the  characteristic  x-ray  blurring 
probability  function,  we  can  calculate  the  number  distribution  of  reabsorbed 
characteristic  x-ray  quanta  blurring  images.  Using  the  characteristic  x-ray  enhancing 
probability  factor,  we  can  calculate  the  number  distribution  of  reabsorbed  characteristic 
x-ray  enhancing  primary  x-ray  input.  With  this  convolution  model,  we  were  able  to 
simulate  the  characteristic  re-absorption  process  inside  a  phosphor  layer  of  an  image 
detector  or  intensifying  screen  by  Monte  Carlo  technique.  The  output  of  the  simulation 
can  be  used  as  two  extra  equivalent  x-ray  inputs,  blurring  characteristic  x-ray  input  and 
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enhancing  characteristic  x-ray  input,  to  every  sub-layer  of  the  phosphor  layer  inside  an 
image  detector.  By  generating  simulated  image  under  an  ideal  point  input,  we  calculated 
the  modulation  transfer  function  of  the  image  detector  with  characteristic  x-ray  re¬ 
absorption.  The  simulation  to  a  Cesium  Iodide  layer  with  pixel  size  17x17  /Mn 2  shows 
that  the  characteristic  x-ray  re-absorption  reduces  the  spatial  resolution  of  the  image 
detector.  The  maximum  reduction  to  the  modulation  transfer  function  of  the  image 
detector  caused  by  characteristic  x-ray  re-absorption  can  as  large  as  10%,  which  is  in  the 
same  range  as  the  measurement.  This  3D  characteristic  re-absorption  model  can  be  used 
to  simulate  the  imaging  process  of  a  phosphor  layer  and  create  various  simulated  x-ray 
images. 
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Figure  1  Schematic  illustration  of  the  generation  and  re-absorption  of  characteristic  x 
rays  inside  the  phosphor  layer  of  an  image  detector  or  intensifying  screen. 


Figure  2  The  characteristic  x-ray  quanta  created  in  voxel  centered  at  (x,  y,  t)  emit 
isotropically. 
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Figure  3  the  relationship  between  the  re-absorption  voxel  centered  at  ( x',y',t ')  of  the  K- 
characteristic  x-rays  and  the  creation  voxel  centered  at(x, y,t) . 
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Domain  B 
(x\y.T) 

Figure  4  The  whole  phosphor  layer  is  conceptually  divided  into  multiple  sub-layers  with 
equal  thickness  At .  The  phosphor  layer  looks  as  a  voxel  set  with  nx  x  N y  x  N,  voxels.  The 

integration  to  calculate  the  total  number  of  reabsorbed  characteristic  x  rays  by  voxel 
centered  at  (x',y',t')  blurring  the  image  is  over  the  whole  voxel  set  except  the  voxels  in 
the  bar  from  the  voxel  centered  at  O',  y',0)  to  the  voxel  centered  at  (x\  y', T) . 
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Figure  5  Schematic  illustration  that  some  of  reabsorbed  characteristic  x-ray  quanta 
enhance  the  primary  x-ray  input. 
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(x',/,0) 
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Figure  6  The  computation  to  calculate  the  number  of  the  reabsorbed  characteristic  x-ray 
quanta  in  a  voxel  inside  the  bar  from  (x\  y',0 )  to  (x\  y',  T)  needs  only  to  cover  the  space  in 
the  window.  The  number  of  total  voxels  involved  is4n2Nr 
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Relotwe  Burring  Kernel  of  the  First  Sub-layer 


Figure  7  Contour  of  characteristic  blurring  probability  function  (table)  on  the  first  sub¬ 
layer  of  a  Cesium  Iodide  phosphor  layer.  The  thickness  of  the  sub-layer  is  2.125  ^ n .  The 
thickness  of  the  whole  phosphor  layer  is  S5/jm.  The  pixel  size  is  17x17  /jin2.  The  index 
should  add  a  bias  -55  when  using  this  2  dimensional  table  to  do  the  convolution.  Every 
sub-layer  has  its  own  table  similar  to  this  table. 
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Figure  8  The  relative  characteristic  blurring  probability  function  of  the  first  sub-layer  in 
a  Cesium  Iodide  phosphor  layer.  A  bias  -55  should  be  added  to  the  index  when  using  this 
2D  table.  Then,  the  center  becomes  (0,  0),  which  is  the  observed  voxel.  The  value  of 
every  element  is  the  relative  probability  that  a  created  characteristic  x-ray  quantum  is 
reabsorbed  by  the  observed  voxel.  The  drawing  shows  that  the  characteristic  x-ray  quanta 
created  in  the  voxels  close  to  the  observed  voxel  have  larger  probability  to  be  reabsorbed 


by  the  observed  central  voxel. 
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Figure  9  Characteristic  x-ray  enhancing  probability  factors  of  all  sub-layers  in  a  Cesium 
Iodide  phosphor  layer.  The  thickness  of  phosphor  layer  is  assumed  as  85  jum  .  The  whole 
phosphor  layer  is  conceptually  divided  into  40  sub-layers  and  each  sub-layer  is  2.125  fjm 


thick. 
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Figure  10  the  number  distribution  of  characteristic  x-ray  quanta  blurring  the  image 
under  the  ideal  point  image  input.  The  distribution  is  obtained  from  Monte  Carlo 
simulation  using  Equation  (20).  The  pixels  close  to  the  central  input  pulse  are  blurred  by 
the  characteristic  re-absorption. 
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Figure  11  the  flow  chart  to  simulate  the  imaging  process  of  a  phosphor  layer  inside  an 
image  detector  or  intensifying  screen 
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MTF  of  The  Simulator  System 


Figure  12  Modulation  transfer  function  (MTF)  calculated  for  a  Cesium  Iodide  layer  from 
simulation.  The  thickness  of  the  Cesium  layer  is  assumed  85  fjm  and  pixel  size  17  x  17//m2 . 
The  phosphor  layer  is  conceptually  divided  into  40  sub-layers.  The  solid  curve  is  the 
modulation  transfer  function  without  considering  characteristic  x-ray  re-absorption.  The 
dashed  curve  is  the  modulation  transfer  function  including  characteristic  x-ray  re¬ 
absorption.  The  spatial  resolution  of  the  detector  is  reduced  by  characteristic  x-ray  re¬ 
absorption 


